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THE LARGE DEFLEXION OF THIN 
CANTILEVERED PLATES 


PART I. GENERAL THEORY 


By J. J. MAHONY (Department of Aeronautics, 
University of Sydney, N.S.W.)+t 


[Received 28 January 1960.—Revise received 30 October 1960] 


SUMMARY 

A scheme of solution, especially suited for digital computers, is presented for the 
case of the bending of a cantilevered plate of small constant thickness, under the 
action of prescribed normal loads which produce moderately large deflexions. It 
takes the form of an asymptotic expansion of the solution of the von Karman 
equations in terms of a large non-dimensional loading parameter. The leading term 
of the expansion is a combination of the solution obtained using the inextensional 
theory of bending and a suitable boundary layer solution. The form of the higher 
order terms is deduced and it is shown how they may be obtained in terms of 
quadratures and solutions of ordinary linear differential equations. 


1. Introduction 
SEVERAL authors have published accounts of an inextensional theory of 


the deflexion of plates under normal load which is well suited for applica- 
tion to thin solid aircraft wings. Mansfield (1) argued that the extensional 
rigidity of a thin cantilevered plate should be so much larger than the 
flexural rigidity that it may be treated as infinite. On this basis, such 
a plate would deflect so that the middle plane of the plate is a developable 
surface. Mansfield (1) and Mansfield and Kleeman (2) showed how the 
application of energy methods enables the developable surface to be 
determined in terms of the solution of an ordinary non-linear differential 
equation. Their work has been applied by Capey (3) and Jones (4) to 
other loading conditions and plate geometries where minor extensions of 
the theory are necessary. Fung and Wittrick (5) observe that if the plate 
were to deflect into a true developable surface it would be necessary for 
a suitable bending moment to be applied to any bent free edge of the plate 
to balance the tendency towards anticlastic curving of the surface. They 
showed how this moment is produced by the action of membrane tensions, 
combined with the curvature of the developable surface, in a narrow 
boundary layer along the free edge in which the plate suffers a small extra 
deflexion from the developable surface. 

Ashwell (6) obtained equations equivalent to those of Mansfield and 

+ Now at Department of Mathematics, University of Queensland. 
(Quart. Journ. Mech. and Applied Math., Vol. XIV, Pt. 3, 1961) 
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Kleeman by considering the equilibrium of a developable surface. This 
work showed that, for general loadings, equilibrium is only possible if 
membrane stresses exist in the plate and concentrated shear forces occur 
on the free edges. The origin of these edge shears is not explained by 
Fung and Wittrick’s boundary layer solution which, however, is not in- 
compatible with their existence. The fact that membrane stresses exist in 
the plate implies that the surface is not exactly developable and in order 
to obtain estimates of the likely accuracy of the inextensional theory 
Ashwell examined an exact solution for the special case of the pure bending 
of a plate of uniform cross-section. From a comparison with the inexten- 
sional theory solution he inferred that, provided the maximum deflexion 
is at least twenty times the thickness of the plate, the latter theory should 
give answers for the deflexion accurate to a few per cent for sufficiently 
compact plates. However, the maximum stress may be underestimated 
appreciably due to the neglect of the boundary layer. Even when the 
maximum deflexion is twenty times the plate thickness the boundary laver 
occupies almost one-fifth of the plate so that it cannot be regarded as 
narrow, an assumption which on the basis of Fung and Wittrick’s treat- 
ment would permit a more accurate calculation of the maximum stress. 
It would appear that the results, obtained on the basis of an inextensional 
surface and a narrow boundary layer, will not be very accurate for the 
ranges of thickness and loading of practical interest in wing problems. 
However, the inextensional theory has great merits in that it is suitable 
for computation on a digital computer and because it can be extended, 
without much difficulty, to plates of variable thickness such as are en- 
countered in wing problems. Thus it was considered worthwhile to seek 
a solution based on the inextensional theory, which is useful for only 
moderately large deflexions and which is not based on any strong assump 
tion as to the width of the boundary layer. As this paper is of an explora- 
tory nature, only the case of plates of uniform thickness is treated in order 
to reduce the algebraic complexity but the methods have been developed 
always with a view to their application to the more general problem. 
In Part I of this paper an asymptotic expansion of the solution of the 
von Karman large deflexion equations is sought, valid for large values of 
a non-dimensional loading parameter. The results of the previously cited 
authors have been combined to give the leading term in this asymptotic 
expansion and then, from an analysis of the properties of this leading term, 
the differential equations and the boundary conditions, the order of 
magnitude of the subsequent terms can be inferred systematically. Substi- 
tution of the formal asymptotic series into the von Karman equations 
leads to two sets of linear partial differential equations. One set of 
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equations describes the higher approximations to the inextensional theory 
and the otherset the higher approximations to the boundary layer. Although 
no real dividing line can be drawn between the boundary layer and the 
developable surface proper under the conditions envisaged here it is con- 
venient to distinguish the two sets of equations as applying to the body 
of the plate and the boundary layer respectively. After suitable trans- 
formations have been made the set of equations for the body of the plate 
can be solved by quadratures, and hence their solution is achieved very 
easily using high-speed computing machines. The boundary layer equa- 
tions are effectively ordinary differential equations with constant coeffi- 
cients and so they also are solvable readily. 

The question of choosing the constants of integration so as to satisfy 
the boundary conditions is not considered in Part I but is given in Part II. 


2. Basic equations 


Consider a plate of uniform thickness ¢ with a typical linear dimension 


a made of isotropic elastic material having a Young’s modulus EZ and 
Poisson’s ratio p. If this plate is subject to a normal load typified by a 
pressure P the equations governing the deflexion of the plate are the von 
Karman large deflexion equations 


VIF = w,- 


u vy? 


femien soe i aie oF w 
Vow = Lp Fig Wyyt Fy Wee — 2 iy Way: 


Wy 


where (ax,ay) are rectangular cartesian position coordinates in the plane 
of the undeflected plate, {12(1—,?)}-*tw(z, y) is the lateral deflexion of the 
plate, {12(1—p*)}-" Z® F(x, y) is the membrane stress function, Pp(x, y) is 
the prescribed load distribution, and the non-dimensional loading para- 
meter L is defined by 
L = {12(1—p*}# Pat/( Et‘). (3) 
For the problems to be considered here, L is a large number and t/a is 
a very small number such that Lt/a is small and the above equations may 
be regarded as having a negligibly small fractional error of order L*#?/a?. 
If s and n are local cartesian coordinates in the plane of the undeflected 
plate which measure distance respectively along the tangent and normal 
to an edge then, on that edge, the boundary conditions to be considered 
here may be expressed in the following form: 
(1) On a free edge 
no bending moment Wants, = 9, (4a) 
no shear forcet Wann t(2—pL)Wyeg = 9,T (4b) 


+ It is assumed that the Kirchhoff expression for the equivalent shear force may be used. 
This implies that the width of the Kirchhoff boundary layer must be small in comparison 


Sass Sg 
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no applied tension F,, = 9, 


no membrane shear force F,, = 0. 
(2) On a clamped edge 
no displacement 


no slope w,, = 0, 


n 


and two of the following conditions: 
no normal displacement Finn t(2+p) Pres = 9, (5c) 
no tangential displacement Fiin—pE i. = 9, (5d) 
no resistance to normal displacement F., = 0, (5e) 
no resistance to tangential displacement F,, = 0. (5f) 


The choice as to which two of the last four equations is to be applied is 
governed by the method of clamping in the plane of the plate. For the 
wing problem conditions (5c) and (5f) are most appropriate because of 
the symmetry. 


3. Form of asymptotic expansion 


There is more than one asymptotic form of solution of equations (1) 
and (2) possible for large values of the loading parameter L. When the 
solution of inextensional theory is applicable the deflexion is proportional 
to the load, which implies that w is O(L) for large values of L. Moreover. 
since the deflected surface is 2pproximately developable the leading term 
satisfies the equation u2,—w,, Wy, = 0 (6) 
and hence equation (1) implies that F is o(Z*). As has been observed by 
Fung and Wittrick (5) and Ashwell (6), the flexural term V‘w alone cannot 
balance the load, except in special cases, so that equation (2) implies that 
F is O(1). It may be noted that equation (6) is a valid approximation 
to equation (1) provided only that F is o(w*) and that the load may be 
balanced if the product of w and F is O(L). Thus there is a range of orders 
of w between L* and L within which the deflected surface will be approxi- 
mately developable. Only the case where w is O(L) appears relevant for 
a cantilevered plate as the work which follows suggests that no mechanism 
exists whereby the relatively large membrane stresses, require 1 for the 
other cases, can be built up with the boundary conditions (4) and (5). 
However, the other limiting forms could occur in problems involving more 
firmly held plates. There is one other possible form in which w is O(L*) 
with that of the boundary layer associated with the developable surface. Later work will 


show that this implies that L+A/a < 1, and hence is satisfied under the assumptions stated 
above. 
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and F is O(L‘) so that the deflected surface is not even approximately 
developable. This form, in which the flexural terms are negligible, applies 
when the plate is clamped on all edges but it may apply also in sub-regions 
of a cantilevered plate under certain circumstances. 

For the remainder of the paper the assumption is made that the defiexion 
is O(L) in the limit as L tends to infinity. Thus it follows that the limiting 
process Jim w(x, y)/L = W(x, y) 
defines a uniformly bounded continuous function W(x, y) which describes 
the developable surface of inextensional theory. It may be noted that 
this is true even if Fung and Wittrick’s boundary layer is included, for 
the additional displacement in the boundary layer is O(1) in the present 
notation. For the membrane stress function the position is slightly 
different because the edge points have to be excluded in the limiting 
process. For all points not on an edge 


lim F(x,y) = F(z, y), 
L->o 


where F,(x,y) is bounded and continuous, but the limit does not yield a 
continuous function at the edge. However, the function W, does not 
provide a uniform approximation to the derivatives of w due to the neglect 
of the boundary layer terms. To overcome this difficulty and obtain a 
suitable leading term in the asymptotic expansion of w it is necessary to 
utilize Fung and Wittrick’s boundary layer solution.t Their result in the 
present notation is that in a region of width O(L~-*) around the free edges 
there is an additional displacement. O(1) and an additional membrane 
stress function O(1). In the limit L > o these additional terms define 
other functions w,(s, L'n) and f,(s, L'n). The combinations 
LW,(x, y)+we(s, L’n), 
Fy(x, y)+-fols, Ln) 

provide uniform approximations to w, F and their derivatives through 
the entire plate as L tends to infinity. For finite but large values of L 
there will be additional terms which will be small in comparison with the 
above. 

To discover the order of these terms it is necessary to investigate the 
boundary conditions and differential equations to see how they may be 
satisfied to successively higher and higher order. Consider the smaller 
order terms which are functions of 2 and y only and whose derivatives are 
of uniform order throughout the plate although not necessarily supplying 


+ For readers unfamiliar with order of magnitude arguments an account of these may 
be found in Goldstein (7). 
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a uniform approximation to w and its derivatives everywhere. If the 
expressions r 
w = LW+A, 


F = F,+B, 


where A is small in comparison with LW, and B is small in comparison 
with F,, are substituted in the differential equations and the effects of 
the boundary layer are ignored, the largest terms in the equation become 


VER, = L{2W,., Ary — Were Ayy— Wey Ace} 


Ory “* ry Orr vy 


and = VA = L{W,,, Byy + Woy Ber—2Wery Bey} + 


Ovy Ory “ry. 


+{A,, FB ytAy% 2A,, F, 


’ 
Orr Ory) 


from which it follows A is O(L-') and B is O(L-*) and a repetition of the 
argument shows that there will be an infinite series in powers of L~? 
However, more terms than the above must be included in order to satisfy 
the boundary conditions and also to allow for the interaction between the 
boundary layer terms and those in the body of the plate. As far as the 
boundary conditions are concerned it has been noted previously that Fung 
and Wittrick’s solution when combined with the developable surface would 
still require an externally applied shear force. An examination of their 
mathematics reveals that this difficulty can be overcome by including a 
further additional deflexion, O(L~') and hence a much smaller one, and 
again confined to a region of width Oj L-') near the free edges. This addi- 
tional deflexion can be caleulated easily by methods closely resembling 
Fung and Wittrick’s but the details are left to section 2 in Part IT. 

The interaction between the boundary layer and the body of the plate 
is estimated most easily by modifying the Mansfield-Kleeman solution to 
allow for the existence of the boundary layer. One of their basic equations 
is obtained by considering the moment equilibrium of any portion of the 
plate on the unsupported side of a generator. For a true developable 
surface only two terms would contribute, the external load and the bend- 
ing moment at the generator. When the boundary layer is allowed for 
there are additional contributions from the membrane tensions in the 
boundary layer due to the additional deflexion from the generator and 
an additional longitudinal bending moment occurring only in the boundary 
layer. It is only the order of magnitude of these terms which is of concern 
in this section so that it suffices to estimate the terms. The two effects 


discussed above will contribute terms whose orders are those of | f,,,,, @, dn 


and { w»,,, dn, and both these are O(L'). Thus, unless these terms cancel 


each other it will be necessary to include a bending moment O( L') in the 
body of the plate to maintain equilibrium. An examination of a particular 
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case, the solution obtained by Fung and Wittrick (5), reveals that the 
cancellation does not occur and there is no reason to expect that this 
result should not be typical. Furthermore, the additional boundary layer 
which supplies the edge shears necessary for equilibrium will likewise lead 
to a term L°W,(x,y). Furthermore, associated with W, and W, there must 
be boundary layers which enable the boundary conditions to be satisfied 
and so it becomes apparent that the asymptotic expansion must be of the 
form a 

> LW, (2,y)+ ¥ L-w,,(s, LAn) (7a) 

) 


m= m= 


and FP = s L-™F (x, y)+ s L-'™f,,(8, L'n). (7b) 
m=0 m=0 

The deduction of the form of this asymptotic expansion is based on two 
limit processes. In one L tends to infinity with x and y fixed. In the other L 
tends to infinity with s and L‘n fixed. The first limit process gives terms 
in the W,, series and the second terms in the w,, series. The process is due 
to Lagerstrom and Cole (8) and has been used by Kaplun and Lagerstrom 
(9, 10,11) and also Proudman and Pearson (12) to solve a mathematically 
analogous problem. The essential difference between the present work 
and that of the cited authors is that the two solutions are superimposed 
instead of regarding the boundary layer and the remainder of the plate 
as distinct regions. However, the work of Lagerstrom and Kaplun (loc. 
cit.) and Proudman and Pearson (loc. cit.) implies that there will be a 
certain arbitrariness, namely the series to which certain terms are to be 
allotted. In this paper the rule is adopted that terms will be allotted to 
that series in which they would appear in the largest order term and when 
doubt remains the terms are allotted to the W,, series. If the behaviour 
of the boundary layer terms w,, are considered for large values of L'n it 
will be seen that any non-vanishing non-negative power of L'n would be 
included in the W,, series by application of the above rule. Thus it follows 
that boundary layer terms will tend to zero for large values of L'n and 
this condition will be shown to determine the division into the two series 
uniquely. 


4. Equations in the body of the plate 


The formal asymptotic expansions (7) may be substituted into the von 
Karman equations and the coefficients of the various powers of L may 
be equated to zero. To each of the equations obtained the limiting process 
with (x,y) kept fixed and not on the edge and L oo may be applied. 
For this type of limiting process the terms in the boundary layer series 
w,,, will not contribute, since w,, and its derivatives, as discussed previously, 
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tend to zero for large values of L'n, i.e. for all points not on the edge in 
the present limit process. The resultant equations are therefore those 
which are obtained by using only the W,, and F,, terms. Thus this process, 
applied to equation (1), yields 


Wee — Woaz Woyy = 9, (8a) 
2WMory Wary Wore Wiyy — Woy Wire = 0, (8b) 
2Wary W, oo” ens ou” Mevy W, p> bn (n > 2), (8c) 


n—1 

, se F a , y 

where dy, mn ae 2, 20h, oy Win— zy — Wo we fa~v yt +n VSF, 4 (9) 
r= 


and ¢, is zero for » less than four and unity otherwise. The result for 
equation (2) is 
Fee Woyy + Foyy Wore —2Fory Wory = VAW— plz, y), (10a) 
Fe Wey t+ Fn yy Wee — 28, cy Macy = Hn (WS 1), (10b) 


nay "Ory 


where n-1 
Pn _ a 2 fa be~ ot gy te oy nti ae 20 ay Wn nay): (11) 


Complicated though these equations may appear they are readily solvable 
once the first approximation W, has been determined using the method 
of Mansfield and Kleeman. . 

The left-hand sides of all these equations, other than the equation to 
the developable surface, are linear expressions in W, or F,, once W, is 
known and, on the right-hand sides, ¢,, involves terms of order (n—1) at 
most while ¥, involves W, and terms of lewer order. The equations are 
therefore suitable for obtaining successive solutions in which first W, and 
then F, are determined for successive values of n. In such a process of 
solution the functions ¢,, and #,, will be known functions at each step and 
so it is apparent that the solution of every one of these equations other 
than for W, hinges on the solution of an equation of the form 


2, Wry —Pee Woyy —Pyy Wore = Gx, y), (12) 


where G@ is a known function. The characteristic lines for this equation 
satisfy the equation 


2 Woy dady + Woy, dy?-+Wo,. dx? = 0 


and it follows from the equation of the developable surface (8 a) that this 
equation for the characteristic directions has equal roots. Equation (12) 
is therefore of parabolic type and the single family of characteristic lines 
satisfy the differential equations 


dy ; » . 
a. 2 rae Woy = —Worz/Woey 
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which, by equation (8a), are equivalent. These equations may be written 
in the differential forms 


Woyy Wy + Wory dz = Woz, dy+M,, dz = 0 


which imply that both W, and MW, are constant aiong the characteristic 
lines of equation (12). However, lines of constant slope on the developable 
surface w = W,(x,y) must be the generators of the developable surface, so 
that the characteristics of equation (12) are the projections on the (z, y)- 
plane of the generators of the developable surface w = W(x, y) given by 
the solution of inextensional theory. The characteristics are therefore 
known straight lines in the plane of the undeflected plate. 

It is convenient to use these characteristic lines to define a new coordi- 
nate system in which equation (12) is reduced to its simplest form. Thus 
let « be any function such that it is constant on any characteristic line 
and such that the value of « serves to distinguish a unique characteristic 
line. A convenient variable a, used by Mansfield and Kleeman, is the 
angle made by the projection of the generator onto the (x, y)-plane with 
a fixed line. In certain cases the value of this angle does not distinguish 
two distinct generators but then another variable is used to calculate the 
inextensional solution and it would appear best to use the same variable. 
The other coordinate, denoted by 8, was determined so that the trans- 
formed version of equation (12) should be as simple as possible for a 
physically satisfactory coordinate system. Extensive algebra shows that 
the optimum choice for 8 is any one in which 8 measures the distance 
along a generator from some fixed point on the generator. Thus here £ 
will be defined as the distance along a generator (or, what is equivalent with 
the accuracy of this paper, along the characteristic) from the edge of 
regression whose tangents sweep out the developable surface. This is a 
convenient origin because the distance f arises quite naturally in the 
Mansfield and Kleeman solution so that the coordinate system is fixed 
very easily once the developable surface w = W,(x, y) has been determined. 
When equation (12) is transformed to this new coordinate system, it 
reduces to the extremely simple form 


Keg = G, (13) 


where A, is the non-vanishing principal curvature of the developable 
surface, and hence is inversely proportional to £8, and so is given by 


Ky wc ko(«)/B, 


where k, is a function known from the inextensional theory solution. 
Thus except on a generator on which k, = 0, ® can be evaluated by 
a double quadrature of known functions, and there will be two arbitrary 
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functions of « which will have to be determined from the boundary 
conditions. The exceptional case of a generator with vanishing curvature 
requires detailed consideration as there is no reason why @ should vanish 
identically on such a generator for all the types of function %,, and ¢,,. 
It would therefore appear that the assumption that the leading approxi- 
mation is a developable surface will lead to a difficulty in this case. This 
is not surprising since the key idea of a developable surface of inexten- 
sional theory is that it is essentially one obtained by a pure bending with 
other deformations being of minor importance. This is hardly a valid 
approximation where the local bending moments vanish. In this paper 
it will be assumed that K, does not vanish. 

It is of interest to note the implications of this solution with regard to 
the form of the deflected plate, apart from the additional displacement due 
to the boundary layer terms. Thus the equation (8b) for W, yields 


W, = hy(x)+-9,(a)B, 


where A, and g, are arbitrary functions to be determined through the 
boundary conditions. It follows that on lines a constant the surface 
w = LW,+ LW, still has a constant slope and hence is a ruled surface. 
For this ruled surface to be a developable surface, other than a cone or 
cylinder, would require that g,(«) = h\(«) since the general solution to 
a developable surface is 


w = h(ax)+h'(a)B. 


However, if one considers the question of determining the function W;, 
following Ashwell’s method and taking the loading as that calculated from 
the boundary layer solution, it is seen that the same number of equili- 
brium conditions have to be satisfied. But the orientation of the generators 
is not now a disposable function, as it was in Ashwell’s case, so that it 
would appear probable that the function g,(«) will need to remain at our 
disposal and not restricted to take the value hj(a). It may be that the 
loads produced by the boundary layer are just right to allow the calcula- 
tion to proceed with one less disposable function but it seems far more 
probable that the second-order approximation is only a ruled surface. If 
this is so then the equation for W, is 


K, Wye x (i), 


so that to this order the surface is curved along the projections of the 
generators of the developable surface of inextensional theory. In any case 
this result will hold once W, is considered unless the membrane stresses 
vanish, which only happens for special loadings. 

In the previous discussion the sets of functions ¢, and y,, have been 





LARGE DEFLEXION OF THIN CANTILEVERED PLATES. I 267 


regarded as known, whereas they are determined from the lower order 
terms according to the operations of equations (9) and (11). However, 
these operations are given there in cartesian form while the solution process 
leads naturally to a solution expressed in characteristic coordinates. Thus 
it is convenient to express equations (9) and (11) in terms of « and f as 
independent variables. For this purpose it is necessary to specify com- 
pletely the coordinate system and here it will be assumed that « is the 
angle made by the generator with tae x-axis. The parametric specification 
of a developable surface implies that both x and y may be expressed in 
the form {g(a)+-g'(a)8} with suitable choice of functions. This together 
with the definitions of « and 8 enables one to determine the transformation 
relations dx = {cos a—f sin a} da+-cos «dB, (14a) 
dy = {sina+-8 cos a}da+-sin adf, (14b) 


and the Jacobian of this transformation has the value (— 8). The case 
when § is zero within the piate is excluded on physical grounds for then 
the external load on the plate must have an infinite moment about the 
generator. When £ is infinite the surface is locally cylindrical and the 
angle « is not a suitable quantity to distinguish generators. In this case 
another characteristic variable is required and the following results are 
not applicable. The transformation relations may be applied to establish, 
for arbitrary functions F and w, that 
rae Cyy TT a Wp, 2Fy Wry 
B-*{ Fy Wee t Fog Wa. — 2F pl Wag +B (wg—w,)]- 
} Fpo{ Bg +B \(wg—w,)|}+-B F,[weg—2w,g—28-\(wg—w,)]- 
+ BF, (1 + B*)wee 2w 8 2B-"weg—w,)| (15a) 


(1+ *)*Feggg + 2B-(1+-28-*— 38-4) Foes 
48-*(1+B-*)Fpgg—B-*(1 —6B-* — 158-*) Fag +-68-*(1 +-38-*) F, 994 
+ 28-*(1+38-*)F ge +B 4(1—6B-*— 158-4) Fz— 108-41 +-38-*) Fg— 
28-9(1+-9B8-*)F, .g—4B*F 448+ 38 (3+ 58) F, + 
+B-*(4+ 15B-*)F,, + 68 Fa t+BPaan: (15b) 


These equations permit the functions ¢, and y,, to be determined readily 
from equations (9) and (11). In practice many of the terms cancel due 
to the properties of W, and W, so that the expression for the early terms 
is much simpler. 

For later work it is convenient to have an expression for the left-hand 
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side of equation (6) and this can be obtained from equation (15a) by 
putting F equal to w. Hence 


Wey — Wes Wy = B*{wag— Wan Wpp—Wppl Prog + B-M(wg—w,)]+ 
+B|2w.g+B-"(wg—w,)|(wg—w,)}, (16) 


and a check on the extensive algebra may be obtained by noting that the 
right-hand side of (16) vanishes when 


w = h(a)+h'(a)B. 


5. Boundary layer equations 


The solution so far obtained neglects the effects in the boundary layers 
apart from the fact that the existence of the boundary layers has been 
used to deduce the form of the expansion applicable to the body of the 
plate. In this section equations are obtained for the boundary layer terms 
in the asymptotic expansion and formal methods for their solution are 
investigated. As the solution in the body of the plate is known more 
conveniently in terms of a and f it is convenient to modify slightly the 
form (7a), (7b) of the asymptotic expansion. This modification is also 
necessary as it is desired to apply the method for the values of J. not very 
large in which case the boundary layers may not be narrow and a coordi- 
nate system based on distance normal to a curved edge may be singular. 
A convenient coordinate system for the boundary layer expansion is («, 7) 
where a is the characteristic parameter as before and 7 is the distance 
along a generator from the edge so that the coordinate system is a non- 
orthogonal curvilinear set. The modified form of expansion to be substi- 
tuted into the von Karman equations is 


w= FLW, («,B)+ & Lm wg(a, 0) (178) 


and F = ¥ LP, (a,8)-+f,(a,9)}, (17b) 
0 


where q a DAs, ; (18) 


8 = Bla)+7, (19) 
and 8 = B(a) is the equation to the edge of the plate and hence B is a 
known function from inextensional theory. In this expansion only one 
boundary layer is included but each generator cuts two free edges so that 
there should be two boundary layer series. Each series may be treated 
independently to a first approximation; the interaction between the two 
boundary layers may be considered afterwards if required. However, 
these terms are likely to be extremely small under most circumstances. 
When expansions (17) are inserted into the von Karman equations, all 
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terms which involve only W,,, F,, or their products cancel out as a result 
of the methods of the previous section. As will be shown, the boundary 
layer terms are exponentially small for reasonably large values of L save 
when » is moderate. Thus if exponentially small terms may be neglected 


B may be replaced by B = Blx)+L-ty 


and the functions W,, and F,, may be expanded in powers of L~* with 
coefficients functions of 7 and a. Thus the result of the substitution of (17) 
into the von Karman equations is a pair of equations involving the sets 
of functions w,,, f,, and known functionst of « and y. If successive powers 
of L-* are equated to zero then two sets of equations result which have 
the form 
A(a) fn ayn = — C(O) m gg t+Om (m = 0,1, 2,...), (20a) 
A(a)0 pq om = C()fnt'¥m (m= 0,1, 2,...), (20b) 
where 
A(a) = (14+ B)?4+4B-%(1+ B)B' +-2B-°(14+-3B-) B+ 
+4B-4*B%4 B-4B", (20c) 
C(a) = B-{Maat+2Mag B’+ B-(14-2B’)(W.+ Mg) + BMg}, (20d) 
and ®,, and ¥,, are functions of the w,, f, with r < m. 
The solution of these equations can be achieved as follows. Consider 
first the equation for the leading terms which may be written 
Sense = —0(a)W,,, (21a) 
Wonnnn = 2(%)fonn> (21b) 
with ao = C(«)/A(a). (21 c) 
These equations are in essence those of Fung and Wittrick’s boundary 
layer and the solution is essentially the same. However, a slightly modified 
argument to obtain the solution is presented here so that the method can 
be applied to the higher order terms. Equations (21) may be integrated 
twice with respect to » and this leads to four arbitrary functions of « as 
constants of integration, all of which must be put equal to zero. If not, 
the solutions for w, and f, would involve cubics in 7 as well as the comple- 
mentary function and such terms will not vanish as » tends to infinity 
as required by the division into w,, and W,, terms as discussed in section 3. 
Thus the integrals of equations (21) are 


Sonn = — a(x), Won aa a(x) fy 
so that Wonynn + OW vais founnn + oP o = 0. 


+ This presupposes that W,, and F,, are known completely without arbitrary functions 
of integration as remain after the work of previous section. In fact, as will become apparent 
in Part II, this assumption does not lead to difficulty. 
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These are ordinary differential equations in » with «a occurring only as 
a parameter and the general solutions are of the form 


exp(yn){C,(a)cos yn +C,(a)sin yn} +exp(—yn){C,(a)cos yn +C,(a)sin yn}, 
where y = {hol}. 


The w,, and f,, terms are required to tend to zero with (x,y) fixed and L 
tending to infinity, that is, » tending to +-co for the case of the edge 
nearer to the edge of regression,+ so that C,(a) and C,(a) must vanish. 
The question of evaluating the remaining arbitrary functions is not con- 
sidered here but is left to Part II]. The same method can be applied to 
solve the equations for w, and f, since the expressions ®, and ‘, are 
exponentially small for large values of » and the argument can be repeated 
ad infinitum. Thus the solution of the boundary layer equations may be 
obtained in principle by the successive solution of a set of ordinary linear 
differential equations. 


+ For the other edge it is 7 tending to — oo. 
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SUMMARY 


The formal scheme of solution, developed in Part I (1), is applied to a particular 
ease for which an analytical solution is obtainable. From an examination of this 
particular solution it is shown how the formal scheme may be completed in the 
general case. Certain limitations on the formal method become apparent and a 
preliminary investigation is made on the possibility of relaxing them. Estimates 
are made of the likely range of usefulness of the asymptotic series in terms of which 
the solution is obtained. 


1. Introduction 

In Part I of this paper (1) a schematic treatment was outlined for obtain- 
ing the solution for the bending of a plate under normal loads when the 
deflected surface is approximately developable. The method was based 
on obtaining asymptotic expansions of the solution of the von Karman 
large deflexion equations for large values of a non-dimensional loading 
parameter. It was shown how the solution could be described in terms of 
two sets of functions each of which satisfied differential equations for 
which general solutions were obtained. The determination of the arbitrary 
functions in these general solutions for given boundary and loading condi- 
tions was not considered and now forms one of the main points of investiga 

tion. One of the difficulties of the method of solution, as far as obtaining 
analytical solutions is concerned, is that the natural coordinate system is 
in general non-orthogonal and curvilinear. This does not affect greatly 
a numerical solution using digital computers, but the associated algebraic 
complexities make it difficult to examine the essential features and limita. 
tions of the method. However, for the special case of a cantilevered 
rectangular plate bent into a portion of a cylinder, the natural coordinate 
system is rectangular cartesian. As this simplification does not produce 
any important change in the mathematics of the method of solution the 
results for this particular problem should be typical of those for the general 

+ Now at Department of Mathematics, University of Queensland. 
[Quart. Journ. Mech. and Applied Math., Vol. XIV, Pt. 3, 1961] 
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case. Any results, which are specific to this particular problem are 
indicated. 

In section 2 the methods of Part I are applied to obtain the first few 
terms of the asymptotic expansion for the deflexion and the stress func- 
tion and methods for satisfying the boundary conditions are discussed. 
It becomes apparent that the form of asymptotic expansion deduced in 
Part I is not complete and the form of additional terms required is 
examined. These additional terms are shown to be significant only when 
the non-dimensional loading parameter is not particularly large. How- 
ever, if these additional terms could be obtained the solution should give 
accurate answers, plus an estimate of the error, for problems where the 
maximum deflexion is only slightly greater than the thickness of the plate. 
A preliminary investigation of possible methods of solution for the extra 
terms is made and their relative importance is assessed. 


2. Rectangular plate bent into a cylinder 

Consider a cantilevered rectangular plate of uniform thickness, acted 
on by a pressure distribution symmetrical about the centre line of the 
plate which is perpendicular to the clamped edge. The notation of section 
2 of Part I will be used with the further specification that the z-axis 
should coincide with the line of symmetry, the y-axis with the clamped 


end, and a, the typical linear dimension, shall be the semi-length of the 
clamped end. It will be assumed that under this symmetrical loading the 
plate deflects symmetrically.t The oniy developable surface which satisfies 
the symmetry condition is a cylinder with its generators parallel to the 
clamped edge. Thus a suitable characteristic coordinate system, as derived 
in Part I, section 4, is the (z,y) system. Furthermore, it is necessary to 
consider only half the plate and this gives some simplification when 
considering the boundary conditions. 

It is convenient to collect here the relevant equations from Part I in the 
modified form applicable to this problem. The differential equations are 


VF = why — Wye Wy (1) (i 1)3 
Ve = Lpl2,y)+ Fez yy +Fyy Wee—2F ry Wey, (2) (Ts 
and the boundary conditions are 

w(0,y) = w,(0,y) = 0 (3a) (I; 5a, 5b) 


+ This is not a trivial assumption. Jones (2) has obtained three solutions, two of which 
are asymmetric for a given symmetrical loading. One can conclude that problems of 
stability such as are discussed by Ashwell (3) are involved. 

t The bracket denotes the number of the corresponding equation in Part I. 
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and two of F(0,y) = F.(0,y) = 0, (3b) (I; 5e, 5f) 


F,,(0, y)—pF,,(0, y) = 9, (3e) (1; 5d) 
F,.2(0,y)+(2+p)F,y(0,y) = 0 (3d) (I; 5e) 

on the clamped end z = 0, 
w,,(x, 1)+-pw,,(2, 1) = 0, (4a) (I; 4a) 
Wyyyl®, 1)+-(2—p)w,,, (2, 1) = 9, (4b) (I: 4b) 
F(z,1) = F(z,1) = 0 (4c) (1; 4e, 4d) 

on the free edge y = I, 

Wrell,y)+pwyy (ly) = 9, (5a) (I; 4a) 
Wrezll, y)+(2—p)we,y ll, y) = 9, (5b) (1; 4b) 
Fil,y) = FAl.y) = 0 (Se) (1; 4e, 4d) 

on the free end x = /, and 
w(x, 0) = w,,, (7,0) = 9, (6a) 
Fi(z,0) = F(z, 0) — 0 (6b) 


on the line of symmetry y = 0. In deriving the equations (3b), (4c), and 
(5c) use has been made of the fact that an arbitrary linear expression in 
x and y may be added to the stress function without affecting the stresses 
derived from it. The form of the asymptotic expansion to be substituted 
into the von Karman equations is 


@ 


0 = 
n 


L1-1°W, (x, y) + s L-w, (x, n), (7a) (I; 17a) 
0 n=0 


x 


F= > L-{F(x,y)+f,(2. 9)} (7b) (1; 17b) 
where pie n = LA(l—y). (8) (1; 18, 19) 
Since the developable surface represented by W, has its generators parallel 
to the y-axis the symmetry condition (6a) implies W, is a function of x 
only. Moreover, since W, represents a surface ruled on the generators of 
W, the symmetry condition implies that W, also represents a cylindrical 
developable surface and is a function of x only. A continuation of this 
argument shows that W, and W, are also functions of x only. As this argu- 
ment depends upon the symmetry of the problem it does not indicate that 
in general the deflected surface will be developable to such accuracy. 

The first few members of the sets of differential equations governing 
the terms in these asymptotic expansions become 


Wo(x)Weyy = 9, (9a) 
We(a)Wy,, = 9, (9b) 


Wo(z)Wayy = 9, (9¢e) 
T 
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Wola) Foy = Wi (x)—p(z,y), (10a) 
Wolt)Fiyy = WY(x)— Wiy(2)Foyy: (10b) 
W(X) Foy = WE (2)— We Boy — Wilz) Fy (10c) 
ete., for the body of the plate and 
fonnan = —W)0Cons (11a) 
Faas ao aS W(2) 19 — Wi(2)Woy n> (11 b) 
Sennan = — Wol®) ann — Wy(2) Cyn — W (2) Conn — Wore Monn + 
+h — Soren (llc) 
Fanny = — Wo(2) a9 — Wy (2) 29 — Wale) 9 — W3(2) Cogn — 
Win More “Onn “122 2fizenm + 2Wory Wiz», etc. ( ll d) 
Wonnnn 42 Wo(®) foun ( 12a) 
iam = Wo(®) fam + Wi(2) foun (12 b) 
Weonnn = W(®) Sonn 7 Wi(2) finn 7 W(X) fonn + Wom Fors zd WorrSonn 7 
+ WonnJorz — 2Woxznn — 2foxn Worn? (12 ¢) 
W3nnnn vit Wol®)fann s Wi(2) fonn > W(2) finn W3(2) fonn wg WorrS run + 
+22 fqn +My Poze +Sorz)— 212299 — 2S ien Worn — 2forn iz 
for the boundary layer. (12d) 
If the asymptotic expansions are substituted in the boundary conditions, 
successive powers of L equated to zero and where applicable the limit 


process L tending to infinity applied to the resulting equations, it can be 
shown that the boundary conditions reduce to 


W,(0,y) = W, (0,y) = 90 all n, (13a) 

w, (0,7) = w,, ,(0,n) = 9 all n, (13 b) 
F.(0,y) = f,(0,n) = 90, all n, (13¢) 

F, (0,y) =f, (0,9) = 0, all n, (13d) 

F, exe(9, y)+(2+p)F, ey, y) = 9, all n, (13e) 


Soran: 0) = frzgn( 9) = (2+B) fn cnn 1)+-Sin—o zzx(9, 9) = 9 (mn > 2) 
(13f) 


and two of 


on the clamped end x = 0, to 
Worn (%, 0) + pWolx) = Wyyp(2,0)+pWy(x) = 0, (14a) 
Wr a lXs 0) + pW, (2) + pen» 2e(%,0) = 0 (mn = 2,3). 14b) 
Wonnn(®, 0) = Wyyyy(Z, 9) = Wongn(®, 0)+(2—H)Wozen(Z, 0) = 0, ete. (14e)F 


+ This result is simpler than the general result since Woz,,(z, 1) vanishes for the solution 
given by inextensional theory as the generator is normal to the edge. Thus the concentrated 
edge shears, introduced by Ashwell (4) do not occur in this case but in general equation 
(14) would require w,,,,,(z, 0) to balance the concentrated shear force. 
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F(x, 1)+f,(2,0) = 90 (all n), (14d) 
fon (2,9) = F, (2, 1) —finsy y(%, 0) = 0 (all n) (14e) 
on the free edge y = 1, and to 
Wil) = 0 (0O<n <3), (15a) 
Won ll, n) es Wryyll, ») = pw, agi’ n) br W(n~2) all, ») = 0 (n a 2) (15b) 
Weil) =0 (0<n <3), (15¢) 
Wornnlls 1) = Wregn (ls 9) = (2—2)Opn cont, 2) + Min-w extn) = 9 (n> 2), 
(15d) 
F(ly) = fill.) = 9 (all n), (15e) 
F, ly) =f, Al.n) = 9 (all n) (15f) 
on the free end x = I. 
The symmetry conditions applied to W, and F,, yield 


n 


W,, (x, 0) = W, yyy(z,0) = 0 (all n), (16a) 


nuyY 


F, (2,9) = Fy yyy{z,9) = 9 (all n), (16 b) 


but a discussion on the conditions on w, and f,, will be left till later. 
Assume now that the functions W,, W,, W,, and W, are known and con- 
sider the solution of the boundary layer equations. Using the arguments 


of Part I, section 5, the equation satisfied by the leading term of the 
boundary layer series is 


—{Wo(x)}*w, (17) 


onmm = 
and hence the general solution of the equation compatible with the pro- 
perties required of boundary layer terms is 


Wy = {Ag(x)cos ky + Bo(x)sin kn}exp(— kn), 
where & is a real positive parameter defined by 
4k4 = {W,(z)}* 
and the corresponding stress function f, is given by 
W (2) fo = 2k*{ — By(x)cos kn + Ag(x)sin kyn}exp(—kn) 


with the same functions A,(z) and B,(x). There are two disposable func- 
tions A,(x) and B,(x) and four boundary conditions, namely, (14a), (14¢), 
(14d), and (14e), to be satisfied. Thus the following conditions must be 


— pWi(x)—2k2B,(x) = 0, 
A,(x)+ By(x) = 0, 

F(x, 1) —2k* B,(x)/Wo(x) = 9, 

A,(xz)+ B(x) = 0. 
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The second and fourth of these conditions are identical and so these 
equations are consistent provided 


F(z, 1) = p (19) 
and the solution for w, is 


Wy = —psign Wo(x){cos ky —sin kyn}exp(— kn). (20a) 
fo = —pi{cos kn +sin kn}exp(—kn). (20 b) 
This is the solution which has been obtained by Fung and Wittrick (5). 


A similar procedure may now be used to determine the solution for w, 
and f,. The integrated form of the differential equation is 
finn = —Wow,— Ww, (21a) 
Wim = Wot, +Wife: (21 b) 
and hence w, satisfies the fourth order linear differential equation 
Wrenn = —{Wo}?1+ Wi fonn— Woo) 22a) 
which, because of equation (11 a) may be written 
’ nie Aas . 7" 7" 4 9. 
Wiongn = —4k*w,—2W, We we. (22b) 
Thus once again there will be only two disposable functions and, as 
above, four boundary conditions to be met. Thus two equations of com- 
patibility must be satisfied and these may be shown to be 


F,(z,1)=0, K(z,1)=0 
by manipulation of the boundary conditions (14a, c, d, e), the differential 
equation (21b), and the properties of the known solution for wy. The 
boundary conditions which determine w, then reduce to equations (14a, c). 


and hence it may be shown that 
w, = pWyexp(—kn){kn(cos kn —sin kn) +-(1+k-*)cos ky 
—(1—2k-*)sin ky}. (23) 
A similar argument may be applied to the solution for w, and after some 


algebra the equations of compatibility of the boundary conditions may 
be shown to be 


0m = 


Ws F(x, 1) LWoxz(X, 0)4 {Wore Sonn + Wonn torr }dndn (24a) 
0 ») 


and Wo F,, (x, l ) MWorry (2; 0) t [ (Worrtonn + WonnJocz} dy. (24 b) 
0 
The process may be repeated systematically for all higher orders to obtain 
conditions on F,(x,1) and F,,(z, 1). 
These conditions may be used to determine the functions W, and F.,. 
Thus if equation (10a) is integrated across the breadth of the plate and 
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use is made of the symmetry condition (16b) and equation (22a) the 
equation 1 


W(x) = | plx.y) dy = Q(x) (25a) 
0 
is obtained, where Q(x) is the non-dimensional load per unit length acting 
on the generator. The result of integrating this equation twice with respect 
to x and using the boundary conditions (15a,c) for the bending moment 
and shear force on the free end is 


z€ L 
Welz) = | | QUE’) de'de = [ QUE'NE’—z) ae’. 25b) 
il z 


The displacement W, is obtained by integrating this « quation twice with 
respect to x and applying the boundary conditions (13a) on the clamped 
end. The stress function F, is obtained by twice integrating equation 
(10a) with respect to y and using equation (19) and either of the equations 
(6b) or (22a). The result obtained is 


uv 


. 
Fy(x, y) = {Wolx)}| wt | | p(x, y') dy’ Qlx)} ay'|. (26) 
1 0 


A similar method can be applied to determine the higher order terms W, 
and F, but the analysis is quite laborious and a simpler method will be 
given later. 

An examination of the above analysis reveals that the only boundary 
conditions on the stress function F, which have been used to determine 
F,, are those on the line of symmetry and the free edge y= 1. The 
boundary conditions (3b, c,d) and (5c) on the clamped and free ends have 
not been used and it is obvious that the solution given by equation (26) 
can at best satisfy these conditions for a particular load distribution. 
Moreover, at the free end, where W(x) vanishes, not only will the boundary 
condition on F not be satisfied, but the in-plane stresses will become infinite. 
Thus the deflected surface will not be approximately developable and the 
whole solution process will fail near the free end. These points will be 
examined further in the next section. 

A method for simplifying the analysis can be obtained by noting the 
physical interpretation of equation (25b). It is the condition of moment 
equilibrium of the portion of the plate to the right of the generator. It is 
important to distinguish between the equilibrium of the plate and the 
equilibrium of a pure developable surface. The latter would yield a 
different condition because allowance for the concentrated edge shears, 
demonstrated by Ashwell (4), would have to be made. The difference 
between the two conditions will be the contribution from the boundary 
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layer and as has been shown in section 3 of Part I this will be O(L*). 
This suggests the simplification of obtaining W,, from the moment equi- 
librium equation of the entire plate rather than using the compatibility 
conditions on the boundary conditions as above. Thus, for example, the 
equation for W, becomest 
Wi(2)+u | wom da— | fom ody = 0 
0 0 
and hence it follows that 


i(%) = —dp(1+4y)sign Wo(x){Wo(z)}'. 


This equation involves calculating the first term w, in the boundary layer 
series, while the previous method involved considering w, and w,. A similar 
advantage persists if higher order terms are sought. A further advantage 
is observed if the previous method is applied to determine W,. The 
equation corresponding to equation (25a) may be shown to be 


W'*(x) = 0 (28a) 


for the case of uniform loading while the moment equilibrium solution 
for this case is 
Wi(x) = —}p(1+4y)sign Wo(x)(l—z). (28 b) 


While this is consistent with equation (28a) it is not the solution which 
would be obtained from that equation by applying the boundary condi- 
tions (15a) and (15c). 

The explanation of this lies in the failure of the solution method near 
the free end as indicated above. This failure is due to the fact that the 
boundary layer on the free edges spreads along the free end as has been 
suggested by Fung and Wittrick (5). This boundary layer can produce 
concentrated tractions on the free end of the developable surface in a 
similar manner to the production of shears on the free edges. Such tractions 
must be considered when applying equation (28 a) which is the equilibrium 
condition for the developable surface. However, they do not have to be 
considered when using equation (28b) for this is derived from the equi- 
librium of the plate (including all boundary layers). It should be noted 
that the two equations are equivalent provided the boundary layer on 
the free end produces an effective shear stress }u(1-+-4,)sign W(x) on the 


developable surface. The moment equilibrium equation is thus seen to 
provide a much simpler determination of the W,, as it does not involve 
considering the boundary layers on the free end. 


+ Strictly the upper terminals of these integrals should be L+. The error in the above 
form will be discussed in the next section. 
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3. Estimation of the error terms 


The above solution has several unsatisfactory features which will be 
considered briefly in this section. It has been observed that not all the 
boundary conditions on the clamped and free ends have been satisfied and 
this indicates the existence of boundary layers on the ends as well as the 
edges. Moreover, Fung and Wittrick (5) have pointed out that the solu- 
tion of the boundary layer equations on the edges is incompatible with 
the conditions on the ends so that the solution must be modified further 
in the corners. There are also errors of the type introduced in the deriva- 
tion of equation (27). Here the aim will be to obtain estimates of all the 
errors introduced by these features and also to indicate possible methods 
for reducing the errors. The arguments used are all of the standard type 
of order of magnitude assessment of boundary layer theory (6). As it is 
only the results of these arguments which are germane to the present 
discussion the detailed argument has been suppressed. 


The boundary layer on the clamped end x = 0 is required to provide 
the in-plane tractions necessary for equilibrium. The application of the 
standard arguments used in boundary layer theory indicates that these 
tractions can be generated in a boundary layer whose thickness is O( L-*) 
and that the leading term in such a boundary layer will be governed by 
the linear equations F 


one = —LW3(0)W,y, 
Worne = LWi(0)F,,. 


Boundary conditions may be obtained from equations (2) and the fact 
that the free edge can produce no in-plane tractions. The solution of 
these expansions can be obtained by expanding in eigenfunctions and 
the expansion is rapidly convergent so that it appears possible to include 
this pattern within the framework of the above solution. The order of 
the additional deflexion in this boundary layer depends on the manner of 
clamping. A leading term O(L~') is necessary in order to satisfy equation 
(3b) while terms O(L~!) and O(L-') are necessary to satisfy equations 
(3c) and (3d) respectively. Away from the clamped edge these additional 
deflexions will be reduced by a factor of the form exp(—const L4z) so that 
the overall contributions of this boundary layer will be small. In addition 
there will be a much narrower boundary layer of width O(L-*) in the 
corner which is necessary to annul the displacement in the boundary layer 
along the free edge. This boundary layer satisfies non-linear equations and 
could be solved by numerical means only. 

The boundary layer on the free end x = / is required to balance that 
portion of the loading which cannot be balanced by the flexural term. 
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It can be shown that this implies that there must be an additional dis- 
placement O(L*) in a boundary layer whose breadth is O(L-*). Both the 
displacement and boundary layer breadth are relatively large so that the 
neglect of these terms is a serious matter. As the equations governing 
this boundary layer are non-linear and would appear soluble only by 
numerical means this term would appear to place the most serious overall 
limitation on the accuracy of the theory. The error introduced by neglect- 
ing such terms will be most serious when the maximum deflexion is being 
calculated as the terms are most significant at the free end. The errors 
are not so serious for the calculation of maximum stress which in general 
will occur near the clamped end, and there the additional term is estimated 
to be Of Lt exp(—const L4/!)]. For general loading conditions the end will 
not be a generator and this difficult additional term will only be required 
to start from the corner of the plate where the deflexion is greatest. In this 
case the above error estimate has to be modified to O[ L* exp(—const L*I)]. 
It sheuld be noted that neither of these error estimates is sensitive to the 
value of L for any but very large values of L but both are much more 
sensitive to the value of /. It is difficult to assess how important these 
terms are without obtaining the values of the constants involved. 

The other form of error introduced, namely, the replacement of L! by 
© in the terminals of integrals of equation (27) and the failure to satisfy 


the symmetry conditions on the boundary layer terms w,, and f,, is seen 
to involve an error which is Ofexp(—const L')]. This is probably a much 
smaller error than the ones considered above but this will depend upon 
the value of / and the appropriate constants. The author has investigated 
the possibility of obtaining such terms and while it is possible it would 
appear that the amount of work involved would not be justified. 


4. Conclusions 


With the error estimates obtained in the previous section it is now 
possible to attempt to assess for any given value of L which terms are 
likely to be important and which may be neglected. The relative im- 
portance of the various terms depends on the aspect ratio of the plate 
but for aspect ratios exceeding unity, i.e. / > 2, it does not seem un- 
reasonable to assume that the exponential terms arising in the boundary 
layer terms may be neglected for values of L exceeding 10. Under these 
circumstances the methods of Part I might be applied reasonably without 
considering the exponential error terms. However, it would not appear 
reasonable to obtain more than three or four terms in any circumstances. 
For values of L near 10 one might anticipate that by retaining the W, 
term in the body of the plate, answers accurate to a few per cent may 
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be obtained. The actual accuracy will depend greatly on the size of the 
coefficients for such values of L. Thus in the example treated in section 2, 
which was the case used by Ashwell (4) to assess probable accuracy of 
inextensional theory, greater accuracy can be obtained with fewer terms. 
This happens because W, is quite small on account of the factor » but this 
is not a typical result for the general problem. An examination of Fung 
and Wittrick’s (5) solution for the boundary layer, when the generators 
are inclined at an angle (47—@) to a free edge, reveals that the solution 
for w, is proportional to (u+-tan*@) and that the thickness is increased by 
a factor sec*#@. Both these effects will increase the size of W, if 6 is much 
different from zero. Thus it is seen that the number of terms in the 
asymptotic series considered in Part I required to obtain a given accuracy 
will vary with the loading pattern. 

Once the exponential terms are negligible, all error estimates are ob- 
tained in terms of the value of L and it is desirable to attempt to correlate 
this parameter with quantities of greater physical significance. A quantity 
used by Ashwell (4) is the ratio y of maximum deflexion to thickness of 
plate which will depend on the loading pattern. To give a definite working 
basis the ratio x has been evaluated for a uniformly loaded plate on the 
basis of inextensional theory and the value obtained is 0-04/*L. It is 
apparent that y is particularly sensitive to the value of / so that the value 
of y is difficult to correlate with the accuracy of inextensional or the 
present theory which is closely associated with L. In fact general criteria 
for the accuracy obtained by stopping the asymptotic series at a given 
stage are difficult to state and are of dubious applicability. There seems 
to be no satisfactory alternative to that of considering each case on its 
own merits. 


For values of L which are not sufficiently small to justify the complete 
neglect of exponentials but for which nevertheless the value of L is not 


too small, the scheme of solution proposed above may still be the most 
efficient. Jones (2) has investigated the alternative procedure of solving 
equations (1) and (2) by an iterative method based on solving equations 


of the form Tv) — 1 


and found that considerable precautions have to be taken to obtain con- 
vergence of the iteration for any value of L which may be considered large. 
His method may be applied to the boundary layer equations for the free 
end and as these do not contain a large parameter convergence difficulties 
should not arise. Moreover, the present scheme involves fewer parameters 
and no parametric functions but, against this, a considerable amount of 
analysis is involved in setting up the scheme proposed here. The decision 
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as to which is the most efficient method for values of Z which render the 
problem difficult would appear to depend upon the number of loading 
patterns and plate geometries to be considered. 
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SUMMARY 

This paper discusses the distribution of stress in the neighbourhood of a pair of 
coplanar Griffith cracks when these cracks are subjected to internal pressure. This 
is a ‘three-part’ boundary value problem and its solution is given by determining 
the coefficients in a pair of dual trigonometrical series. In the special case of uniform 
internal pressure, closed formulae are given for the resulting shape of the cracks 
and for the stresses in their neighbourhood. These formulae involve the well- 
tabulated elliptic integrals and some numerical results are given. 


1. Introduction 

Tue stress distribution in the interior of an infinite two-dimensional 
elastic medium when the very thin internal crack —c < y <c (x = 0) is 
opened by internal pressure has been studied by Sneddon and Elliott (1). 
Their solution depends on that of a pair of dual integral equations. When 
the medium contains the two coplanar Griffiths cracks —b < y < —a 
(zg = 0) and a <y <b (x = 0), Sneddon’s expressions for the stresses 
and components of the displacement vector lead to a set of triple integral 
equations. By making use of the Weber-Schafheitlin integral, these in- 
tegral equations can be solved by determining the coefficients in a pair 
of dual cosine series. Formulae for finding these coefficients have been 
given recently (2), and the ‘two crack’ problem is theoretically solved in 
this way. 

Detailed consideration is given to the special case in which the internal 
pressure is constant. In this case, closed expressions are obtained for the 
normal stress and displacement in the plane of the cracks and numerical 
results are given for various values of the spacing ratio a/b. 

The method of the present paper can be considered as a direct extension 
to ‘three-part’ boundary value problems of that devised by Sneddon and 
Elliott for ‘two-part’ problems. The solution given is an alternative to 
that proposed by Willmore (3), who based his work on the complex 
variable methods developed by Green and Stevenson and by Muskhelish- 
vili (4). 


(Quart. Journ. Mech. and Applied Math., Vol. XIV, Pt. 3, 1961] 
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2. The triple integral equations 

The stress in the medium can be described by three stress components 
Fz, Gy, Try If the displacement vector is (u,v, 0), the boundary conditions 
are that wu, v and all three stress components vanish as x? y? tends to 
infinity, and that +, = 0, o, = —p(y), the internal pressure applied to 
the cracks, when zx = 0, —b <y <a anda <y <b. Symmetry about 
the plane z = 0 shows that, in fact, T,y vanishes for all values of y when 
x = 0, and that u = 0 on that part of x = 0 which is unoccupied by the 
cracks. For simplicity, p(y) is taken to be an even function of y so that 
there is symmetry also about the plane y = 0. 

Sneddon (5) has shown that, in z > 0, suitable expressions for o,, and 
u are 


2 o 
ee } x(p)(1+-ap)e~*? cos yp dp, (1) 
0 


ie aL { x(p){2(1—o) -+ap}e-ze SP dp, 
7 p 
0 


where Y, o are respectively Young’s modulust and Poisson’s ratio for the 
medium and x(p) is a function to be found from the boundary condition 
on « = 0. Sneddon also gives expressions for o,, ,,, and v, but these will 
not be required in this note. 

When z = 0, equations (1) and (2) give 


x 
. 


9 
[or]e-0 = —= | x(p)eos yp dp, 
WJ 


0 


4(1—o%) f 
=> = | x(p) 2H? dp, (4) 


[w],-9 = yet aces 
0 


and the boundary condition on x = © leads to the triple integral equations 


x 


x(p) dp = U (0 i y < a) 
0 P 

| x(p)cosyp dp = {xp(y) (a<y <b) 
0 


@ 


| x19) SM? ap (y > b) 
: 


+ Young’s modulus is denoted by Y instead of the more usual Z to avoid confusion 
with the elliptic integrals which appear in the solution. 
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3. The determination of y(p) 
If n is an integer, particular cases of the Weber—Schafheitlin integral 
are (6) 
cos{(2n— 1)sin-'(y/b)} 
rs 2__)i 
} Jy,,-,(bp)cos yp dp = im si 
(y?—b*)Hy + (yb) 


x | cos{(2n— 1)sin *(y/b)} 


(y < 6), 





(6) 





(y > b), 
"i “ene EE < b), 
| Jon s(bp)—? dp = n—I » ates (7) 
: 0 (y > b). 
Hence by taking x(p) ein y A, de (bp). 
n=1 


where A,, are coefficients to be determined, the third of equations (5) is 
automatically satisfied. This choice of y(p) also ensures that the value 
of |u|,» is continuous at y = 6, a condition that is clearly required on 
physical grounds. 


Writing y = bsin }@, a = bsin $f, (9) 


the first two of equations (5) become, after interchange of the order of 
integration and summation and use of (6), (7), (8). 


vA 


> An 5, 008 3(2n -1lé=0 (0<6< B) 


n=1 


>A , cos 4(2n—1)6 kb p(bsin $0)cos 36 (B<0<7) 
Taking 

B.=6 fie 40 p(bsin $@)cos }(2n—1)0d0  (n = 
so that : 


> B,, cos 4(2n—1)0 = $b p(bsin }0)cos}@ (0 < 6 < 7-7), 
n=1 
equations (10) can be written 


— cus }(2n —1)0 = f(@) (O<0< B) 


s (A, —B,)cos}(2n—1)0=0 (B<O0< nm) 
nel 


x 


where f(0) = — > m. 7 cos 4(2n —1)0. 


n=l 


NEC MACS CE AVON Se Ceremer 
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The coefficients in the dual cosine series (13) are given by (2) 
4(A,,—B, )eosec 38 = &(1)F(n, 1—n; 1; sin? $8)— 


1 
— | &'(8)F(n,1—n; 1; 8? sin? 48)ds, (15) 


0 





du +-constant, (16) 


AG ie ae 
where &(s) = ~ | uf" s “nee $8)} 
0 


the constant being determined from continuity ia the value of [w]|,_, at 
y = a. It is also worth noting that the value of > (A,,— B,,)cos }(2n—1)@ 
n=l 


when 0 < 6 < f is given by (2) 
1 
&(1) [ £’(s) ds 


abi aie coionstomelontia ee — 17 
(1—sin? $6 cosec? $8)' a (s*—sin? $6 cosec? $8)! mu 
sin $8 cosec } 





Once the coefficients A,, have been determined from equations (11), (14), 
(15), and (16), the function y(p) is found from equation (8). All the stress 
and displacement components can then be found from the expressions 
given by Sneddon and the problem is theoretically solved. 


4. The coefficients A,, in the special case p(y) = p, (constant) 


For the special case in which the cracks are subjected to uniform 
pressure, p(y) = pp» (constant) and equation (11) gives 


7 


B,, = bp, { cos 4.cos }(2n—1)0 dd = io 3 : “ (18) 
0 at a: 


Equation (14) reduces to f(#) = —4abp, cos $0, so that f’(@) = }abp,sin 40 
and f’{2sin-(usin 48)} = }rbp,pysin}f8. Hence, from (16), 


x 2 j 
[ a + constant = }bp,(s?+C)sin}f, (19) 
6 — pL . 


&(8) = bp,sin $f 


» 


) 


C being independent of s. 


Writing ssin}8 = sin }A and expressing the hypergeometric functions 
as Legendre polynomials (7,a), equation (15) yields 


A,—B P 
_ = (1+ C)sin? SBP,,_,(cos B)—} | »-y(cOSA)sinAdA. (20) 


0 





ON COPLANAR GRIFFITH CRACKS 287 
5. The determination of the constant C 


The constant C in equation (20) is determined from the continuity of 
[w},.9 at y= a. From equations (4), (7), (8), and (9) 


a =a le u],-0 “>< in ey 008 (20 — 10 (y < 5). (21) 


n=1 


Using (18) and (20), after interchanging of the order of integration and 
summation, this can be written 


_Y[u),-0 | ‘ 


Osi — P,_,(cos B) 
ss +4(1+ C)sin? ps batt. on—1)0— 
4bp9(1—o*) Be nadia Zi a Ghee niece ad 


B . 
oo P,,_,(€08 A) ‘ 96 
“a4 sin A i E008 (20 — 18 dd. (22) 
0 = 


The sum of the series in equation (22) is found in the Appendix to be given 
F(cos 4A, }7) (0 < A), 


{ (ore +} . 
“1 4A, sin cost A) | (@ > A), 


Pe ee A) 


cos $(2n—1)@ = 


where. in the usual notation for elliptic integrals, 


sind 


F(k,¢) = far dd i =| seit Gl ae 


~k? sin?¢)! f(1—#?)(1—k*2?)}!3 


It is also shown in the Appendix that 
F Jcos 42, sin (e een 


—) on 7 ai ‘cos $0 21 | (cia) 
cos 40— E jos $8, sin~ {= )} 1B) + sin SBF cos $8, sin- os 8)! 


] 


é sod 
wheve E(k,¢) = | (1—k?sin2$)! db = [ (Se Je dt. 
a Sete 
Continuity in displacement requires that [w],_, = 0 when y = a, that 


is, when 6 = 8. For these values of [w],_, and @, equations (22), (23), and 
(25) give 


0 = }(14+ C)sin? 48 F (cos $8, 47)+4E(cos $8, 47)—4 sin? 48 F (cos $8, 47), 


leading to 3(1+-C)sin? 48 = sin* 1b eee. (27) 
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6. The displacement when a < y < 6 


The shape to which the cracks open under uniform internal pressure 
Py is given by equation (22) with 8B < @ <7. Writing 
cos 3 b?—y*\4 . 
28 
mee cos $8 =(5=5) ; (8) 


and using equations (23), (25), and (27) we find 


Fisl.. E(cos $8, 47) 5, 4 
poll —ot) ~ E(cos 48, 6) — Nios oth te 5 F (cos 48, 4). (29) 


It is interesting to note that when a = 0, the cracks merge into a single 
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crack of width 2b. In this case, equation (9) shows that 8 = 0 and 
cos$8 = 1. The elliptic integrals degenerate into 
¢ y?\4 
E(1,4¢) = | cosd dd = sing = ud 
0 
d 
aa sec ¢ dd = 


and, since F(1,¢) > 0 as ‘ > 47, equation (29) reduces to 


gp ltee, aoe (1 a ai 
2bpy(1—o? b? 
This is the formula given by Sneddon for the single crack and it shows, 
of course, that the crack opens into an ellipse. 
Fig. 1 shows the shapes of the cracks for various values of the spacing 
ratio a/b (= sin $8). These curves are easily calculated from equations (28), 
(29) and tables of the elliptic integrals E(cos 48, 4), F(cos 48,4). As was 
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to be expected, the cracks become ‘egg shaped’ when a/b is small and 
approach ellipses as a/b increases. 


7. The normal stress when 0 < y < a 
It also seems of interest to obtain a formula giving the normal stress 
co, in the plane of the cracks. From equations (3), (6), (8), and (9) we find 


_jnb[o,],-9 = sec 40 4, cos }(2n—1)8 


= 4nbp,+sec 40 > (A, —B,,)cos}(2n—1)0, (30) 
n=1 
the last step resulting from the use of (18). Now the range 0 < y <a 


12 
ror 
0-8 
0-6F 
0-4 


0-2} 








| ehh 
0 02 04 O06 O8 10 
Fic. 2 





corresponds to the range 0 < @ < B and use of equations (17), (19) gives, 
for such values of 0, 


B,,)cos 4(2n—1)0 


1+C 28 ds 


—sin? $0 cosec? $f)! — . {¢—sint 40 cosec? $8)*| ” 
sin $6 cosec 4 
(31) 


Substituting for (1+-C) from (27) and evaluating the integral on the right, 
the right-hand side of equation (31) reduces to 


E (cos $8, $7) 


F(cos $8, 4) 
‘ 








twbp, sin iB), ] 


}nbpo|sin? 40— }(sin® 3—sin* 40)-1. 
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Substituting in equation (30) and using (9), we have 
lozle-9 _ 5*{ E(cos $8, $7)/F(cos $8, $r)}—y* 
Po {(a®—y*)(b?—y*)}# 
The stress at the origin is given by 
pf Blcca df. 47) sal 
P*"\aF (cos 48, $m) |” 


and its variation with the spacing ratio a/b is shown in Fig. 2 





—il (®9<y<a 
(o, |, 0-0 = 


8. The normal stress when y > 6 


When y > 6, the stress [o,],_, is given by equations (3), (8), and (6). Use 
of equations (18), (20), (27) and the result (7,5), namely 


> (—eafts (fay =3(%j—sintp)', ( 
oa, 1) ea 1) P,,_,(cos A) Wes sin? 4A} , (34) 


gives, after some spaseae 
[oz|-9 _ y*—b*%{ E(cos $8, =) F(cos $8, $7)} ' an 
rat | - b), 3: 
Po {(y*—a*)(y?@—b*)}! Ps Nat 
and this should be compared with equation (32). 
When a= 0, B= 0, E(1,4n) = 1, 


F(1,47)—+ © and equation (35) 
reduces to 


oe pe 
loz|,-0 = Pals 2} (y > 5); 


this is Sneddon’s formula for a single crack of width 2b. 


APPENDIX 


1, We require the sum of the series 


S > Pe 16) co. 4(2n—1)0. 


2n 
n=l 


Starting from the result (7,5) 


> h?*-*P,,_,(cosd) (1 — 2h? cosA+h'*)-4, 
n=l 


integration with respect to h between 0 and z gives 


22/(1 + 2%) F 
1) l P t 
P, ,(cos A) ila eerie — x — A — EE SE s 
l J (1— 2h? cosA+hA*)t fe {(1 —t*#) 1 —¢ cos* 4A)}? 
n=l 0 0 
(38) 
when we write ¢ 2h/(1+h*). Writing z = e*” and equating real parts 
sec 4 
dt 
{(1—@)(1 —#? cos* $A)}* 
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When @ < A, this gives 
' d 
t , 

| aay Reo ~ Pleo ae 

ny 
in the notation of equation (24). When @ > A, equation (39) becomes 
L sec 10 
dt dt re 


oN 7 TS AEE ON olileanaee 2 Ae Oe ‘Saiele 
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writing ¢cos }A = 1/t’ in the second integral. This can be written 
cos £4 see 5A 


it’ 
2S = 2F(cos }A, }7)- - 


) aro py 


0 


. os $6) | 
F \cos $A, sin (— ‘)}. (41) 


using the periodic property of the elliptic integral. Equations (40) and (41) give the 
complete result (23) used in the text. 


2. We also require the value of the integral 
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writing k cos $A. Now, from (26), 
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k k(k? — cos* 40)! 
the argument of the elliptic integrals having been omitted for brevity. From 


equation (24) 
(1/k)cos 40 
| dt 


[oo | er 
F\k, sin '{ 7, cos 48), “1 ey - kt?) yh" 
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giving, (7,¢), 
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dk (—e@fhl—eey 
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where k’? = 1—k*. This reduces to 
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when use is made of equations (43) and (44). Hence (42) gives 
I = 4[{E—(1—*) Frag: 
Using (26) and (24) 
cos 46 
E{1, sin-(cos $6)} } dt = cos $0, 
¢ 


03 48 4 , 19 
t + COs 
T'S) aint ‘ = ee stpionapeaicresinaaot 
F{1, sin-(cos $9)} | ‘eer } log, ——— ia) 
0 
and hence (1—k*)F vanishes when k = 1. Using these values at the upper limit in 


equation (45), we find the value for J given in equation (25). 
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SUMMARY 

It has been shown by Synge (1) that the propagation of Rayleigh waves on a free 
plane surface of a transversely isotropic medium is only possible if the free plane is 
either normal or parallel to the direction about which there is symmetry of rotation. 

In this paper, the displacements due to surface waves radiating from a given 
harmonic source are found in terms of double Fourier integrals. In order to satisfy 
the boundary conditions, additional free vibrations of the semi-infinite medium have 
to be introduced, These are expressed as double Fourier integrals, the integrands of 
which contain poles and branch points. 

Using the isotropic case as an analogy, the branch points of the integrands are 
neglected and the double integrals are estimated asymptotically, using a method 
derived from the one given by Lighthill (2) for triple integrals. 

With this technique it is shown that in the first of the two cases suggested by Synge 
the slowness and wave curves are circles, and that the asymptotic forms of the 
displacements are, consequently, of a simple nature. The resulting equation for 
the velocity of Rayleigh waves is the same as that obtained by Stoneley (3) and other 
writers by the method of plane waves. 

The second case is fully anisotropic, when the method of plane waves is insufficient. 
Equations representing the slowness and wave curves are obtained. These curves 
are plotted from data obtained numerically for a number of materials. Calculations 
show that in certain circumstances, in this case as well as in the first case, the variation 
of wave amplitudes with depth from the free surface is harmonic, in addition to being 


exponentially damped. 


1. Introduction 

THe geometry of the propagation of three-dimensional waves in aniso- 
tropic media is best understood in terms of two surfaces, the slowness 
(or phase) surface and the wave surface. The two surfaces are polar recipro- 
cals with respect to a given sphere and can be explained physically in the 
following way. If O is the centre of the given sphere, S is a point on the 
slowness surface, and P the corresponding point of the wave surface, then 
the vector OS represents the direction and reciprocal of the magnitude of 
the velocity of waves observed at P emanating from a given source at O. 
Contributions to the wave motion at P are made by all points S; at which 
the outward normal is parallel to the OP direction. The wave surface 


itself represents a surface of equal phase of waves at time ¢ = 1 which 
Now at Dept. of Applied Mathematics, The University of Sydney, N.S.W., Australia. 
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started from O at time t= 0. If we make use of ideas suggested by 
Lighthill (2), the vector OP represents the group velocity of energy propaga- 
tion from the source in the OP-direction, while the corresponding vectors 
OS, represent the inverse phase velocities of the waves which are propa- 
gated in the direction OP. In general, OSP is not a straight line, except 
in the isotropic case, when all the surfaces reduce to concentric spheres. 

In considering infinite, anisotropic, elastic media, Musgrave (4) used 
harmonic analysis into plane waves to obtain the slowness surface whence 
the wave surface was obtained by a method which considered the energy 
flux across surfaces of equal phase. After expressing displacements in 
terms of Fourier integrals depending on the radiation from a given harmonic 
point source, Buchwald (5) used a new method developed by Lighthill (2) 
to estimate asymptotically the multi-dimensional Fourier integrals, and 
obtained all the geometrical properties of the waves, as well as asymptotic 
expressions for the wave amplitudes, in terms of the distance from the 
source and certain invariants of theslownesssurface. The connexion between 
parabolic points of the slowness surface and cuspidal edges or conical 
points of the wave surface was also explained, and it was shown that the rate 
of decay of waves in directions corresponding to these singular points of the 
wave surface is less than in other directions. Duff (6), using another Fourier 
integral method, has made a detailed analysis of the Cauchy initial value 
problem and showed that there are sharp waves (corresponding to the 
wave surface), radiating from a given source, together with continuous 
waves of smaller amplitudes which are usually confined to the space 
between the sharp waves corresponding to the outer and inner sheets of 
the wave surface. 

A full account of the geometry of Rayleigh surface waves in anisotropic 
media can only be given in terms of one or other of the above-mentioned 
methods. At the cost of some extra labour, a Fourier integral method is 
used in this paper so that information concerning the amplitudes of the sharp 
waves in addition to the geometry of their propagation can be obtained. 
Whilst only the special case of transverse isotropy is considered, the methods 
used can be employed for surface waves in all anisotropic media, and need 
not necessarily be confined to applications in the field of elasticity. 

Synge (1) considered Rayleigh waves in anisotropic elastic media as 
the limiting case of an excited surface layer. He arrived at a number of 
general conclusions, including the fact that, except in special circum- 
stances, undamped Rayleigh waves do not exist in anisotropic media. It 
can be shown that one condition for the existence of undamped Rayleigh 
waves is that the free plane surface corresponds to a plane of symmetry 
of the anisotropic medium. 
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In this paper the simplest anisotropic case is considered, namely trans- 
verse isotropy, where there is rotational symmetry perpendicular to a 
given direction. In such a medium the above condition for the existence 
of Rayleigh waves is satisfied only if the free plane is either perpendicular 
or parallel to the given direction, normal to which there is symmetry. Both 
these cases are considered in this paper, although only the latter is truly 
anisotropic. In the former case all directions parallel to the free surface are 
equivalent so that as far as Rayleigh waves are concerned the medium is 
quasi-isotropic. This case will be used as a first simple example of the 
method, since results may be compared with the equation for the velocity 
of Rayleigh waves which was obtained by Stoneley (3) and other writers. 
The solution of the more complicated and fully anisotropic second case is 
dealt with in the latter part of this paper. 


2. The radiation from a source in a transversely isotropic medium 
If (u,v, w) is the displacement vector referred to a rectangular system of 
cartesian coordinates: O(x,y,z), then the equations of motion are 
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where p is the density and X, Y, Z the body force resolutes. 
In a homogeneous, transversely isotropic material the stresses can be 
expressed linearly in terms of the displacements by the stress—strain 
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when it is assumed that displacements are sufficiently small for Hooke’s 
law to be satisfied. There are five independent elastic moduli, and it can 
be shown that the stress-strain relations, together with the equations of 
equilibrium, form a system which is axially symmetric about the z-axis, 
corresponding physically to isotropy at right angles to the z-direction. 
In order to take advantage of the axial symmetry we express the dis- 
placements in terms of potential functions defined by 

ad dp 
+ Bt 4 
ey eax” 
20 


oz 


@ oz 


and 


It is to be noted that ¢ rolb nn ge ss as so that % corresponds to the 
ex by ex y* 

z-component of the rotation. Eliminating the stresses and displacements 

from the above equations we obtain, after some reduction, the three 

equations of motion 


29 23 r a 

e é oY ox 
a,—._.+a, Vi——/V? - —— = 

ez? et ox cy 


A2 2 
ne nee tes : oi = 0, 
ez 


et2}az? * ez 
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2 

and a,Vi-4+ (“65 =4+a4V ag Vv? 


where Vj is the two-dimensional Laplacian operator é*/éx*+-é?/éy?, and 
the five constantsa, = ¢,,/p,@_ = Cgg/p, 2g = (Cag +Cy3)/p,%q = (Cy: —Cy2)/2p, 
Gs = Cy/p, have the dimensions of (velocity)?. It is important to note 
from (2.6) that % is independent of the other displacements, whilst ¢ and @ 
are coupled by the equations (2.7) and (2.8). 

We assume that a harmonic source of circular frequency w can be 
represented by body forces of the form 

(X,¥,Z)=e- [ [ [ (X,Y, Zexpfi(ar+py+yz)] dadBdy, (2.9) 
so that the Fourier transforms 


(X,¥,Z) = — =o Lid (X,Y, Z)exp[ —i(ax+ Py+yz)| dadydz 
© —o (2.10) 
are known functions of a, B, y. 
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Let /, 4, @ be expressed as the integrals 
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d = e-tat | | | dexp|i(axr+By+ yz)| dadBdy, (2.12) 
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and A=eto [ [ [ Dexpli(ax+fy+y2)] dadBdy. (2.13) 


—o— rs 


Substituting the expressions (2.11) and (2.9) in (2.6) we find that the fune- 
tion w& in (2.11) is given by 

b = fy/Fa,B,y), (2.14) 
where F(a, B, y) = a,(a?+-8?)+a,; y?—w?, (2.15) 
and (a?+-f*) fy = i(BX—aY). (2.16) 
Similarly, from (2.7), (2.8), (2.9), (2.12), and (2.13) it can be shown that 


db = fs/G(a,B,y), 2.17) 


where 


G(a,B,y) = [a5(a?+ 8?) +a, y?—w?][a,(a? + 6?) +4, y?—w*?|—aj y*(a? +B), 
2.18) 


(a?-+ 8%) fy = tay y(a2+f*)Z —ifa,(a?+6*)+a,y*—w"(aX+BY), (2.19) 
and that 8 = fe/G(a, B, y), (2.20) 
where y*fo = tay y*(aX +BY) —iy[a,(a?+f*)+4, y?—w? |Z. (2.21) 

For example, an isolated harmonic force acting in the z-direction at the 
origin can be represented by 

X=Y=0, Z= Z,d(x)d(y)d(z)e-™, (2.22 


where 4(2) is the Dirac delta function. In this case it can easily be shown 
that 


fg = fg = iagyZ,/8n°; and fg = —ija,(a*+f*)+a,7*—w?]Z,/8n*y. 
(2.23) 

3. The quasi-isotropic case 
Here we assume a semi-infinite medium in the region z > 0, the plane 
0 being a traction free boundary. The next step is to reduce by one 
dimension the integrals in (2.11), (2.12), and (2.13) using a contour in the 
complex y-plane consisting of the real axis together with a large semicircle 
in the upper half-plane, the contributions from this semicircle being vanish- 
ingly small (by Jordan’s lemma). With this device, the integration with 
respect to y can be replaced by appropriate contributions from the singu- 
larities of the integrands. Since the functions f,, f,, and fg represent a 
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source at the origin they do not have any singularities except, possibly, at 
infinity and at either « = 8 = 0, or y = 0. However, these latter poles 
are due to the introduction of potential functions for the displacements 
and it can be shown that any contributions to the potentials from the 
singularities at « = 8 = 0 and y = 0 yield zero displacement everywhere. 
Thus the only important singularities of the integrands are at F = 0 in 
(2.11), and at G = 0 in (2.12) and (2.13), where the functions F and G 
are defined by (2.15) and (2.18) respectively. Given real values of « and f. 
@ is a second order polynomial in y’, whilst F is linear in y*. The real 
roots of F = 0 and @ = 0 have been discussed in a previous paper (5) 
and correspond to ordinary waves propagated in three dimensions, the 
single sheet of the real surface F = 0 and the two sheets of the real surface 
G = 0 making up the three sheets of the slowness surface of these waves. 
For Rayleigh waves we require that the amplitudes should decay exponen- 
tially with distance from the free surface. It is easily seen that this condi- 
tion is satisfied if we look for roots of F = 0 and G = 0 for real values of 
a and 8 such that no real positive values of y” satisfy either of these equa- 
tions. This is so when a, 8 are made large enough in order that the line 
with these values of «, 8, parallel to the y-axis, does not intersect either 
of the real surfaces F = 0 or G = 0. 

With these conditions satisfied, the equation F = 0 has a pair of conjugate 
roots for y on the imaginary axis, whilst the equation G = 0 has two pairs 
of conjugate complex roots which all lie on the imaginary axis if the 
discriminant of the quadratic equation G = 0 for y* is positive. If this 
discriminant is negative the four roots are complex, each — a real 
part and an imaginary part of the same magnitude. Let y = +ir be the 
two roots of F = 0, r being real and positive, and y = +-is,, +78, be the 
four roots of G = 0, where s, and s, are either both real and positive, or 
are complex conjugates with positive real part. Replacing the integration 
with respect to y in (2.11), (2.12), and (2.13) by the contributions from the 
appropriate poles included in the semicircular contour defined above, 
the radiation from the source in the form of surface waves is given by 


ys = Qrie— to! { ae exp[t(ax+ By)—rz| dadp, (3.1) 


ee Y 
$ = 2nie~ot 5 Ey Ei a exp[i(ax+y)—s,2]dadg, (3.2) 
hes a one - 


2 9 


and @= dmie-iet J Je | exp[i(ax+By)—s;z]dadB. (3.3) 
=1,3 yuts; 
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Here F, = 0F /éy, G, = @G/éy, and r, 8,, 8, are the roots of 
a,r? = a,(a?+f?)—w’, (3.4) 
and a, a, 8*— Hs*+-[a,(a?+-B*)—w?*|[a,(a?+-B?)—w*] = 0, (3.5) 
respectively, where 
H = (a, a, +a? —a§)(a?+f*)—(a,+a,)w?. (3.6) 
The boundary conditions on the face z= 0 are zz = yz = 22 = 0. 


From (2.4) and (2.5) it can be shown that the three conditions on #, d, and 6 
at z 0 are op 


= 0. (3.7) 
: (p+80) = 0, (3.8) 
ez 


v2 6*6 3.9 
and (a,—as) 1o+as—; = 0. (3.9) 


Clearly the functions given by (3.1), (3.2), (3.3) do not satisfy the boundary 
conditions, and we must postulate additional free vibrations 


x 


b* = 2rie~tot ( [ A(x, B)expfi(aa+By)—r2] dadp, (3.10) 


x 


zx x 


mie [ [ Bia,B)expli(ae+fy)—s;z]dadB, (3.11) 


and * = 2mie-to#t ( [ C,(a, B)exp[i(ax-+By)—s,; 2] dadB, (3.12) 
where j takes the values 1, 2. The additional potentials ~*, ¢*, and 6 
satisfy the equations (2.6), (2.7), (2.8) with zero-body force if r, s; are given 
by (3.4) and (3.5), and 

a;(a*+-B*) B;+[a5(a?+B?)—a,s?—w?®|C; = 0 (j = 1,2), (3.13) 
where the five functions A, B,, C; are such that the functions %+-~*, 
}+-¢$+¢43, and 6+-67+-6} satisfy the boundary conditions. 

From (3.1), (3.7), and (3.10), ~+%* = 0 everywhere in order that the 
boundary conditions may be fulfilled. We conclude that the function % 
makes no contribution to Rayleigh waves in this case. The functions 
6+¢}+ 43 and 6+6$+6} must satisfy (3.8) and (3.9), whence 


§,(B,+C,)+s,(B,+C,) = Pr, (3.14) 
ys «| 45" ~ (3.15) 
Stab Fy Synte 
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Writing 6? for «*?+-f? and a, for a,— 
a, 6*( B,+ B,)—a,(s} C,+83C,) = A, (3.16) 
where A+ a aed —a [2 | (3.17) 
2 ; y= tay 1G yates! 


j~=1,2 ¥ 


Equations (3.13), (3.14), and (3.16) form a set of four equations in the 
B,; and C;, whence 


B, = win (3.18) 
and C; = (3.19) 
where the functions g,(«, 8), A,(«, 8) are related to the source terms I’, A by 
hy(a, B)/ag = b*g,(a, B)/(a,6*—a, s3—w) 
= Tb*{a,(a, b?@—w*)+a, a, 83|— As,(a,6?+-a, 83+), (3.20) 
with similar equations for g, and h,. The function K(b*) is given by 
K(b*) = 8,[a,6?+a, s{-+-w*|[a4(a, b?—w*) +a, 4; 83] — 
~89| ag b? +-a, 83 -+-w? |[a,(a, b?—w*)+-a, 4, 87]. (3.21) 
Taking out the factor s,—s, and using the sum and products of the roots of 
(3.5), this expression reduces to 


K(b*) = (8,—8,){(a,b?—w*)|[ (a, a, —ag)b*—a, w?|—a,a,8,8,}. (3.22) 


In the particular case of the isolated harmonic force acting in the z- 
direction at the origin, [ and A can be shown to be given by 


ae 2) a, 8? +a, b? —w* (3.23) 
"823 8,(4a, rr er 2H) see 


=1,2 
and = 7 PS (545-04, )0" —Gy(y 8} +e? ) (3.24) 
8;(4a, a, 87 —2H) 


In this case [ and A are functions of 6?, but for a source which does not 
possess rotational symmetry about the z-direction, [ and A would not be 
axially symmetric and g,; and h; would not be functions of 5’. 

Inspection of (3.5) shows that G = 0 and G, = 0 cannot hold simul- 
taneously when y = is,, whence it may be seen that, apart from the point 
at infinity, the only singulerities of g; and h; are branch points at s, = 0 
and s, = 0, and a possible pole at b = 0, which may be neglected. 

We may conclude that the important singularities of the integrands in 
the expressions for ¢* and 67 are poles at K(b*) = 0, and branch points 
at s, = 0 and at s, = 0. 

Before carrying out the next step of estimating asymptotically the 
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double integrals for the potential functions, using the method described 
in the appendix, we should note two important points. 

(1) If the expressions (3.2) and (3.3) are examined, it can be seen that 
the integrands do not contain any important singularities (since G = 0 
and G, = 0 cannot hold simultaneously) and, therefore, the primary 
waves from the source make no significant contributions to the waves 
at infinity. 

(2) In using the semicircular contour mentioned in the appendix, the 
branch points s, = 0 and s, = 0 are included within the contour taken, 
so that this contour will have to be indented along cuts to exclude these 
points. The integration along the cuts gives extra contributions to the 
integrals. However, in considering isotropic media it has been shown 
by a number of writers (7) that while the poles of integral expressions 
obtained correspond to Rayleigh waves, the integrations along the cuts 
represent other waves due to the free surface which decay as O(R-*), 
compared with O(R-*) for Rayleigh waves. 

Taking the isotropic case as an analogy, the contributions to the asymp- 
totic estimation of the integrals for ¢7 and @f from poles corresponding to 
K(b*) = 0 represent Rayleigh waves, whilst the contributions due to the 
existence of the branch points represent other kinds of waves. 

We are now left with the following four expressions: 


* g(a, B) io 
| KO) exp[i(ax+By)—s,;z| dadB (3.25) 


¢; = 2wie~* 
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(Ms K(b?) P) expfi(ax+By) —8;2| dadB, (3.26) 
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" 2 

and OF — Qie—iot { 
x 


where j = 1, 2, and the double integrals are to be estimated asymptotically 
for large values of R = ,/(x?+-y*) by the method described in the appendix. 
We must remember, however, that the integrations along the branch cuts 
are being neglected. 

Using the nomenclature adopted in the appendix the slowness curve is 
given by K(b?) = 0. Taking out the factor s,—s, from (3.22) the curve is 
given by 

(a; b?—w*){(a, a, —a§)b* —a, w*}—a, a, 8, 8, = 0. 
Using the product of the roots of equation (3.5), we find, after some algebra, 
that 


a, M*b*’— M(M +-2a,a,)b‘w* + (2M +a, a,—a,a,)a, b*w*— 
—d,(d,—a,;)w* = 0, (3.27) 
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where M = 4a,(a,—a,). If we write c = w/b, this equation is identical to 
the one obtained by Stoneley (3), Sato (8), and Das Gupta (9) and, when 
a, = a, = ci, a, = ch, ag = a,—a,, it reduces to the standard equation 


ee ee 


for the velocity of Rayleigh waves in isotropic media. Stoneley showed 
that the cubic in (3.27) has only one real root, which is less than va;. In 
consequence, K(b*) can be written in the form 

K(b*) = 1(b2)(c6*—w), 
where L(6?) is a quadratic expression with no real factors, and w/c, the 
positive root of equation (3.27). 

Clearly the slowness curve in this case is a circle of radius w/c, and the 
wave curve a concentric circle of radius c,. The fact that these curves are 
circles is, of course, due to the quasi-isotropic natare of this case. The 
curvature of the slowness curve is cp/w, and |grad K| = wL(w?/c},)/cp, 80 
that the expression |grad K| |\k,|* in expression (A.8) of the appendix 
reduces to (w/c,)'L(w*/c},). There is, of course, only one point on the 
slowness circle where the normal is parallel to a given direction of propaga- 


tion, and sgnk, = —1, so that the asymptotic forms of the integrals for 
d; are 


; ~ (2m) pee Eat xp[i(apx+Byy—wt+ }r)—s;z], (3.29) 


with similar expressions for 6;, the only difference being that h,(«, By) is 
substituted for g,(a»,8,). The point (ao, 8) is on the slowness circle, so that 
od + 6% = w*/c}. The constant L(w*/c3,) is, in fact, taken care of by equal 
constant factors in the numerators g(a», By) or h,(a», 8), which appear on 
solving the equations (3.13), (3.14), and (3.15). 

If we now substitute b = w/c, in the equation (3.5) for s, we can find the 
decay factors s;. The equation (3.5) has either two roots which are real 
and positive, or two complex roots which have equal positive real parts 
and imaginary parts which are equal in magnitude but have opposite sign. 

The results of calculations on a number of materials are given in Tables | 
and 2, in all cases the magnitudes of elastic constants are those given by 
Hearmon (10). In Table 1, the first column gives the velocity cp, of Ray- 
leigh waves, and the second column the ratio cp/,/Ja;, which corresponds 
tO Cp/¢, in isotropic media. For the materials listed in this table, the roots 

8, of equation (3.5) are real, and in columns 3 and 4 the decay coefficients 
‘ aaa o, are given, where oc, = 27¢,,8,;/w, and the decay with distance from 
the free surface is given by the factor e~-%*, ) being the wavelength. 
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In Table 2 are listed two materials, zinc and apatite, for which equation 
(3.5) has complex roots. Columns 1 and 2 are as for Table 1, but columns 
3 and 4 now give op, a, the real and imaginary parts of c, so that the decay 
corresponding to both roots of (3.5) is the same but the phases have opposite 
sign, the decay factor having the form exp[—(e,+i;)z/A|. It is most 
interesting to note that periodic decay can occur for media of this type as 
compared with purely exponential decay in isotropic media. Zinc and 
apatite are also the most interesting cases encountered in section 4, where it 
will be shown that for these materials the decay is purely exponential, or 
periodic as well, depending on the direction of propagation. 

Table 3 gives the values of the five elastic constants, calculated from 
numerical data given by Hearmon (10). 


TABLE 1 





Material CR 


/ 
RINGS 





Cobalt 
Magnesium 
Beryl 

Ice 
Beryllium 


2°80 km/sec 
2°91 Pm 
47° 

1°78 


j 7°25 
i 





0-962 
0°940 
0956 
0958 
0926 














TABLE 2 





Material 


€R 


ER/NGs 





Apatite 
Zinc 


3°94 km/sec 
2°06 





0858 
0°868 











TABLE 3 
Elastic constants (after Hearmon, 1956) 





Vaterial 


ay 


a, a; 


ay 


as 





Cobalt 
Magnesium 
seryl 

Ice 
Beryllium 
Apatite 

Zinc 








47024 
3454 
8-752 
1°629 
19°53 
4°446 
o’704 


2°080 














| in units of 


| 10"! cm?/sec? 





4. The anisotropic case 

In this case the free plane surface is parallel to the direction about which 
there is rotational symmetry (the z-direction). Without loss of generality 
we assume that the free surface is the plane y = 0 and that the medium 
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oceupies the space y > 0. The boundary conditions are now 
yy =ty=2y=0 aty=0. 
From (2.4) and (2.5), these reduce to 
x (+0)—% = (4.1) 
(2. ay 6 
dy? azx* + éxdy 
e o0 Op 
oy? :|$+ = io exey 


(4.2) 


ge 
la. 2a,)- at an (4.3) 


where a, = a,—a;. 

Using a method similar to that used in section 3, it can be shown that 
the radiation from a harmonic point source in the form of waves on the 
surface y = 0 is given by 


2 a 


} = Qrie JJ ] exp[i(axr+yz)—ry| dady, (4.4) 
Ff, ir 


== Qrie—tt ps j | exp[i(ar+yz)—s;y]dady, (4.5) 
j=1,2 _*, "B41 B= tay 


and 6 = 2mie- > | J aR expli(ar+yz)—s,y|dady. (4.6) 
a is; 


j=1 


These expressions are obtained by performing the integration with respect 
o B in the triple integrals for y, 4, @ in (2.11), (2.12), and (2.13) where, as 
in the previous section, the integration with respect to 8 is replaced by the 
residues at the poles at F = 0, G = 0 respectively; F and G being given 
by (2.15) and (2.18). The numerators f,, f,, and fgare given by (2.14), (2.17), 
(2.20) respectively and Fz; = ¢F/@8, Gg = @G/aB. The quantities r and s, 
are here the imaginary or complex roots of the equations F = 0, G == 0, 

for 8, respectively, so that r is the positive real root of 
a,r? = a,a*+a,;y?—w?, (4.7) 


and 8,, 8, are the two real and positive (or complex with positive real part) 
roots of 


a, a,(8*— «*)*— H(s?—a*)+ (a, y?—w?)(a; y?—w?) = 0, (4.8) 
in which H = a,(a, y?—w*)+a,(a, y?—w?)—a3y*. It is to be noted that 
if the roots of (4.8) are complex, then all four roots have equal real parts 


and equal imaginary parts, and we select s, and s, to be the two complex 
conjugate roots whose real parts are positive. 
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The above expressions for the primary radiation from the source do not 
satisfy the boundary conditions, and we must add ad ditional free vibrations 


@ «© 


b* = 2rie~™ [ [ A(a,y)expli(ax+yz)—ry] dady, (4.9) 
$* = 2mie~ [ [ Bia, y)expli(ax+yz)—syy]dady, (4.10) 


and 0? = 2aie~ = ( ( O,(a, y)exp[i(ax + yz)—s,y] dady, (4.11) 


where r and s, satisfy equations (4.7) and (4.8) so that the functions ~*, ¢F, 
@* satisry the equations of motion (2.6), (2.7), and (2.8) and, therefore, 
from the last of these equations 
a,y°C;+ [a5 y*+4,(a?—s7)—w? |B, = 0. (4.12) 
The functions J+*, 6+¢f+¢3, 6+0f+0F must satisfy the boundary 
conditions (4.1), (4.2), and (4.3), whence 
8,(B,+C,)+8,(B,+C,)+ tad of (4.13) 
(r?+-a*?)A —2ia(s, B, +8, B,) = A, (4.14) 
and > (a, — 2a,)a® —a, 87|B;+ 2ia, aA —agy*(C,+C,) = A. (4.15) 


j 


The functions [', A, A are given in terms of the source functions by 


r+ ial fs AA > ed 0. (4.16) 
? 1 Bir ae Gg \p-is 


=1,2 


A+(r?4 ha 4 a - (0). (4.17) 
Fg }\g-ir Pm a] Btw 


til. 


: Se | Le 
I (a,—2a,)o2 a] A gy" = 
a ‘ Gg B is; Gg B iss) 


Equations (4.12), (4.13), (4.14), and (4.15) form a system of five simul- 
taneous linear equations for the five unknown A, B; and C;. Solving, we find 


and 


A(a,y) = Ray (4.19) 


By(a, y) = a in B. (4.20) 


' hj(a, y) 
and O(a, y) = ren (4.21) 


002 55 xX 
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The numerators p, g; ,and h; can be found in terms of I’, A, A, and the 
coefficients of the five equations (4.12)-(4.15), whilst K is given by 


rat+oa® Qos, 2as, 
K(a,y)=| dgay* J, 82], |, 
2a,a,ar Jd, J | 
where, in this expression, 


I, = agy*? —a,(a?—87)+w? 
and J; = 2a, a3 0*—a, a;(a?—87)+a,(a, y?—w?). 


By arguments similar to those used in section 3, the primary functions 
Ws, 6, and @ make no contributions to the waves at infinity, and the important 
contributors to these waves are the singularities of the functions A, B,, C; 
which, from (4.19), (4.20), and (4.21) are branch points at r = 0, s, = 0, 
and s, = 0, and possible poles at K = 0. Assuming that any contributions 
due to the branch points correspond to waves which are not Rayleigh waves, 
the slowness curve for Rayleigh waves is given by the equation 


K(a,y) = 0, (4.23) 


in which r is the positive root of (4.7) and s,, 8, the appropriate roots of 
(4.8). This equation is homogeneous in a, y, and w, so that the phase 
velocity w/,/(«?-+-y") is independent of w and there is no dispersion. A close 
inspection shows that (s,—s,) and (a, y*—w*) are factors of the function 
K(a,y), but that the poles associated with these factors give trivial contri- 
butions to the displacements. The determinant in (4.22) could, therefore, 
be simplified by the removal of these two factors. However, the system of 
equations is very complicated and since we have to resort to numerical 
solutions in any case, there seems little point in performing this minor 
simplification. One important point must, however, be made at this stage. 
The quantities s,, s, are either both real and positive, or complex conjugates 
with positive real part. It can be shown quite simply that even in the latter 
case, removal of the factor s,—s, reduces K = 0 to an equation with real 
coefficients since only terms containing 8, 8, and 8,-++-s, are then involved. 

Ignoring any contributions due to the existence of the branch points, 
the asymptotic forms of the potentials are, using the technique given in 
the appendix, 


aay 2 Ye) 
b~ as > : A, a ae exp|i(a, x t+Y¥q2z—wt) ry|, (4.24) 
ai"q 


q 1 


with similar expressions for ¢; and @;, for ordinary points of the slowness 
curve K(a,y) = 0, and with expressions derived from that of (A. 12) of 
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the appendix for points with zero curvature. From (A. 10) of the appendix 
the wave curve of Rayleigh waves is given by 


z 1 
8 ee (4.25) 
uK, +yK, 


y 
in terms of the parameters a, y which are connected by the equation 
K(a, y) 0. 

The curves given by (4.23) and (4.25) are too complicated to discuss 
analytically, and it was necessary to resort to numerical methods to plot 
the slowness and wave curves for particular materials. The slowness 
curves obtained are given in Figs. l(a), 2(a), 3(a), and 4(a) in which 
the heavier line represents the slowness curve of Rayleigh waves and, for 
comparison, the finer lines represent the intersections of the three dimen- 
sional slowness surfaces in an infinite medium with the plane B = 0. 
Figs. 1 (b), 2(6), 3(6), and 4(6) represent the corresponding wave curves; 
here again the heavier line represents the Rayleigh waves, and the finer 
lines represent the intersection of the wave surfaces in infinite media with 
the plane y = 0. It should be remarked that although we mention only 
the wave curve of Rayleigh waves, in actual fact Rayleigh waves should 
be represented by cylinders parallel to the y-axis intersecting the plane 

y = 0 on the given curve. 

_ From the calculations, the materials considered have, as in section 3, 
grouped themselves into the two sections illustrated in Tables 1 and 2. 
Figs. 1(a) and 1(6) refer to cobalt, and both ice and magnesium have 
very similar properties. The behaviour of beryllium is rather like that of 
beryl, which is illustrated by Figs. 2(a) and 2(6). The main feature here 
is that in a narrow wedge in the x-direction no propagation of Rayleigh 
waves occurs. 

The behaviour of zine and apatite is more interesting. In Figs. 3(a@) and 
4(a) the broken lines represent the value of y* above which equation (4.8) 
has complex roots. Thus, in Figs. 3(b) and 4(b) Rayleigh waves in the 
wedges containing the x-axis decay exponentially in the y-direction in the 
usual way but in those containing the z-axis the decay is harmonic as well 
as exponential. Another feature common to both zine and apatite is the 
occurrence of two cusps of the wave surface close together near the z- 
direction, so that there are thin wedges in which N = 3 in the expression 
(4.24) and three Rayleigh waves are propagated. 

As is shown in Fig. 3 (a), there are no Rayleigh waves in zine with phase 
velocities in the x-direction, although propagation in this direction does 
take place. From Fig. 3(b) it can be seen that this is due to waves whose 
group velocity is in the x-direction, but whose phase velocity is not. 
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(c) 
Fic. 4 


Fic. 1. (a) Slowness, and (6) wave curves of cobalt. The heavy lines represent 

the Rayleigh waves and the finer lines are the intersections of the slowness 

and wave surfaces of three-dimensional waves with the planes y = 0 and 
y = 0 respectively. 


Fic. 2. (a) Slowness, and (6) wave curves of beryl. The Rayleigh waves are 
shown by the heavier lines. 


Fic. 3. (a) Slowness, and (b) wave curves of zine. The thicker curves repre- 

sent the Rayleigh waves. The broken lines in (a) give the values of y* above 

which the variation with the depth y is periodic as well as exponential. In (6) 

the decay with the depth of y of these waves is purely exponential in the 

wedges which are bounded by the broken lines and contain the x-axis, while 
in the wedges containing the z-axis this variation is periodic as well. 


Fic. 4. (a) Slowness, (b) and (c) wave curves of apatite, (c) being a magnifica- 
tion of part of (6). The broken lines have the same meaning as in Fig. 3. 
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Finally, propagation in the x-direction in apatite is very complicated. 
Fig. 4 (c) isa magnification of the relevant part of Fig. 4 (6). It demonstrates 
the fact that several cusps of the wave curve exist close together, and that 
as many as five Rayleigh waves can be propagated in some directions. 
However, it is improbable that five waves as close together as these can be 
distinguished in practice, and all that one would therefore expect to 
observe in this material is a concentration of Rayleigh waves travelling in 
the x-direction. 

An interesting point which arises out of Figs. 3 (a) and 4 (a) is the apparent 
connexion between conical propagation of the three-dimensional waves 
and the existence of a region in the ay-plane in which the equation (4.8) 
has complex roots. This connexion can be confirmed theoretically, but a 
brief description of what is meant by conical propagation is first necessary. 
For example, in Fig. 3(a) the slowness surface of waves propagated in 
three dimensions in zine is obtained by rotating the figure (omitting Ray- 
leigh waves) about the y-axis. The surface so obtained is given by the 
equation G(a,8,y) = 0, where the function G(a, 8, y) is defined in (2.18). 
The planes generated by the rotation of the broken straight lines in the 
figure touch the surface on two circles (called circles of tangency), the 
planes of which are parallel to the plane y = 0. On the wave surface, a 
single point on the z-axis corresponds to the whole of each circle of tangency. 
Such a point is called a conical point, and the corresponding wave propaga- 
tion decays as O( R-*), compared with O( R-') for ordinary waves in three 
dimensions [{(2), section 5 and (5), section 5}. 

Taking G(a,0,y) = 0 as the equation of the intersection of the slowness 
surface with the plane 8 = 0, the parametric equation of the intersection 
of the wave surface with the plane y = 0 is 


ce Y y i tate oe: ’ 
t~= —G, G., z= G, ¢.., 


) 


(4.26) 
At a conical point, z = 0 in (4.26) for some non-zero value of a. Since 
G(a, 0, y) is a quadratic function of a?, the condition for the existence of a 
conical point is that the equations 

ey GoM 

~ Oe) 


Y 


should hold simultaneously for real values of «, y when 8 = 0. In other 
words, for some real y the quadratic equation G(«, 0, y) = 0 in a® must have 
equal roots. The equations in (4.27) can now be solved simultaneously, 
giving the intersections of the circles of tangency with the plane 8 = 0. 
However, if we replace a? by a?—s? in the equation G(a,0,y) = 0, we 


obtain the equation (4.8) for s?. Therefore, when y? is such that (4.27) is 
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satisfied then, at the same time, the quadratic equation (4.8) for (a? —s*) 
has equal roots. When »* = w*/a,;, the equation (4.8) must have real 
roots, and it can also be shown that at any circle of tangency y* > w?*/a;. 
The conclusion is that if a circle of tangency exists, then the value of y* 
giving this circle is the same as the value of y* below which the decay of 
Rayleigh waves with depth is purely exponential, and above which the 
variation with depth is harmonically as well as exponentially damped. 
In Figs. 3(a) and 4(a) the broken lines, touching the circles of tangency, 
are the divisions between the two kinds of decay. 


5. Conclusion 

It is clear that the method of plane waves alone is insufficient to describe 
the properties of Rayleigh waves in anisotropic media. In order to do this 
a method should be used which gives the geometry of actual wave propaga- 
tion and, preferably, some additional information about wave amplitudes 
at large distances from the source. The only case when analysis into plane 
waves yields the correct result is the one considered in section 3 in this 
paper. This quasi-isotropic case has been considered by Stoneley (3), 
Sato (8), and Das Gupta (9), but it should be stressed that their approach 
of considering plane waves in particular directions will not yield enough 
information to describe the geometry of propagation in more general cases. 

There remains the question of whether Rayleigh waves in anisotropic 
media are at all possible. Synge (1) suggested that, in general, they cannot 
be propagated, and showed that in the case of transverse isotropy, Rayleigh 
waves can then only be propagated in the two cases discussed in sections 
3 and 4 above. He pointed out that for arbitrary cases the function K 
which defines the slowness curve is a complex function of two variables. 
The equation K = 0, therefore, is satisfied only at the intersections of two 
curves which meet at isolated points, so that the propagation of Rayleigh 
waves can only take place in directions corresponding to these intersec- 
tions. However, a deeper examination of this problem gives rise to two 
hypotheses. (1) The complex equation K = 0 can have solutions corre- 
sponding to propagation in all directions on the plane surface, but the result- 
ing waves are dissipative except in the directions of the above-mentioned 
intersections, when there is no dissipation. (2) The equation K = 0 is a 
real equation if and only if the free plane surface is a plane of symmetry 


of the anisotropic material. It is hoped to test these hypotheses in due 
course. 
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APPENDIX 
THE ASYMPTOTIC ESTIMATIONS OF CERTAIN 
DOUBLE INTEGRALS 


1 estimating asymptotically the value of integrals of the form 


antes IE PP expliloa By)) dad, (A. 1) 


ox 


Pa fs 
for large R (x? +y?), we make use of a method derived from one Lighthill (2) 
developed for triple integrals. Similar techniques have been used in optical and 
electromagnetic wave problems by a number of authors, quoted by Born and Wolf 
(11), and by Bunkin (12). However, the method is little known and for that reason 
it is worth while discussing in this appendix. We assume that f (a, 8) has no singulari- 
ties except at infinity, and that the only singularities of the integrand are poles 
at K 0. We rotate to new coordinates 2’, y’, where the y’-direction is parallel to 
the vector R whose components are (x, ¥). The above integral reduces to 


x 


e7iwt ( 1 f(a’, B’) 


Ka’. By Py) daa. (A, 2) 


We now carry out the integration with respect to 8’, making use of a contour in the 


complex 8’-plane consisting of the real axis together with a semicircle of large radius 


-— «x 


in the upper half plane. Applying Jordan’s lemma we find that the only contributions 
to I are from the poles of the integrand at K = 0. We can ignore complex values of 
x’, B’, which satisfy this equation, as these values give contributions to J which are 
exponentially small for large R. We are left with contributions from the real curve 
K(a’, 8’) = 0 in the a’, 8’ plane, corresponding for fixed a’ to poles on the real axis 
of the complex f’-plane. It can further be shown that contributions from poles which 
do not satisfy the radiation condition are ignored. The integral in (A. 2) now reduces to 
I ~ Qrie~*+t sé oe xp[ip’y’} da’ [K, Bh (A. 3) 
where a’, 8’ are evaluated on that part of the curve K(a’,f’) = 0 for which the 
radiation condition is satisfied. If we assume that K is a homogeneous function of 
x, 8, and w, then the radiation condition reduces simply to B’ > 0 in the special 
coordinates, and to ax+By > 0 in the original coordinates before rotation. 

Making use of the method of stationary phase, the chief contributions to the 
integral in (A. 3) are from the neighbourhood of those points on the curve K(a’,B’) = 0, 
where 8’ > 0 and where the normal is parallel to the B’-direction. Let there be N 
such points, with coordinates (a,,8;), at which 8’ can be expanded in the Taylor 


series B’ = B+ bk, (a’—a;,)?+..., (A. 4) 


where k, is the curvature. 
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4) 


Using the result f exp(fiku®) du — J 


7 )exldarisen ki 
we find that pw S 4 —S1%Bp) ripe y’—wt)) 

Yt La Klay Bibel 
where the phase constant A, is given by 


A, = exp[}ri(2+sgnk,)]. (A.7) 


The integral in (A. 6) can now be expressed in terms of the original coordinates by 


f (ap, By) : 

ty ad KL, [et Plat + By yo] (A.8) 
where |grad K| = ./(K2+ Kj). The curvature k, can be expressed in terms of the 
derivatives of K (indicated by suffixes) as follows: 


a K2 Kgg— 2K, Kg K,,+ Kp KK, 


ky (K2+ K3)t 


(A. 9) 
The sign in the expression (A. 8) is positive if the function K is increasing and negative 
if K is decreasing with distance from the origin at the point (a,,f,). 

If the integral J in (A. 1) represents waves from a source situated at the origin, 
then the waves observed at a distant point, represented by the position vector R in 
(x, y) space, correspond to those points on the curve K(a,8) = 0 where the normal is 
parallel to R. 

The curve K(a,8) = 0 is called the slowness or phase curve, and the contributing 
points on this curve give the directions and inverse phase velocities of waves at the 
point of observation. Clearly R — OVK. 


where © is a scalar, so that if a represents the vector whose coordinates are (a, 8),then 


R = ——.. VK. 
a.VK 
If we give the term wave curve to curves of equal phase at t = 1 which started from the 
origin at t = 0, a.R = w, so that 
w0K VK 
et es a am rere (A. 10) 
a.VK eK /ew 

is the equation of the wave curve in terms of the parameters a, 8 which are connected 
by the equation K(a,8) = 0. In this equation we have assumed that 


a.VK —w0eK /éw. 
This is true only if K is a homogeneous function of a, B, and w, when there is no 
dispersion. The slowness and the wave curve are polar reciprocals with respect to a 
given circle whose centre is the origin. 


The case of zero curvature of the slowness curve must also be considered. Let 
k, = 0, so that the Taylor series for f’ is 


B’ . B, r th, (a’— a)? + see (A. l 1) 


Using this expression in the integral in (A.3) we find the contribution from such a 





RAYLEIGH WAVES IN TRANSVERSELY ISOTROPIC MEDIA 317 


point to the integral is given by 


fr 4 r 4 ‘ a 
I, ~ ry) (4) Oe Po explilfy 9 —wt)]. (A. 12) 
Similar expressions can be found if the first non-zero term of the Taylor expansion 
for 8’ is of higher order. 
It is to be neted that waves corresponding to a point of inflexion of the slowness 
curve decay only as O( R-+), compared with O(R-4) for ordinary points, and it can 
be shown that there is then a cusp on the wave curve. 
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SUMMARY 


Some general and particular solutions of the electromagnetic equations are given 
for steady flows with finite electrical conductivity. In most cases the electromagnetic 
forces are assumed to be too small to perturb the flow significantly, but some of the 
solutions are also dynamically exact. Both plane potential flow of an incompressible 
fluid and rectilinear flow are considered, and the magnetic fields are imposed by 


current-carrying wires. The equations are derived for both plane and axisymmetric 
conditions. 


1. Introduction 


We suppose that a stationary magnetic field is provided by insulated 
current-carrying wires within a fluid of finite electrical conductivity. If 
the velocity field and the variation of electrical conductivity are regarded 


as known functions, the magnetic field can be determined from Maxwell's 
equations and Ohm's law alone; following Lundquist (1), we term this the 
kinematic problem. When the magnetic force per unit area is negligible 
compared to the stagnation pressure of the fluid, the velocity and magnetic 
fields can be determined separately as a fluid dynamic and a kinematic 
problem respectively, When the magnetic forces are small but not 
negligible, we may usually treat the first approximation for the magnetic 
field as a kinematic problem. 

Solutions for magnetic fields in free space can be derived from the vector 
potential A of the magnetic field given by 


Aw tte 


J P 
where I is a line current vector, dl the length of a small element on the 
line current, and p the distance from that element to the point where A 
is evaluated. We shall show that analogous methods can be used for 
certain types of kinematic problem. As in classical theory, idealized con- 
cepts, such as line and sheet currents, are used. Physically a line current 
would correspond to a fine wire whose dimensions are too small to perturb 
the flow significantly. A sheet current might be given by a series of closely 
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spaced parallel wires, and in general it will be assumed fully permeable 
to the flow. It is convenient mathematically to represent the applied 
currents, which are in practice carried by the wires, as a general current 
density field j, assumed to be known. We suppose that the medium which 
carries j, is so tenuous that its physical presence can be neglected, and, in 
spite of any applied current field, a flow region is completely filled by fluid. 
j, is allowed to exist artificially at any point in the flow, and is independent 
of j, the current density within the fluid. From an analysis based on this 
idealization we shall be able later to derive equations in which the applied 
currents take a physically more acceptable form, i.e. line or sheet currents. 
The magnetic field depends on the total current density at a point 


curlH = j+jo, (1) 

where the m.k.s. system of units is used. In a steady state 
curlE = 0. (2) 
Also div B = 0, (3) 
B = pH, (4) 
and the permeability » is assumed to be constant. Within the fluid Ohm’s 
law gives j = o(E+v~ B), (5) 


where o is the electrical conductivity; we shall assume that o is constant 
throughout the flow. 

Kinematic problems can be reduced to their simplest form when the 
following conditions are satisfied within the fluid: 


Case A(i). In cylindrical polar coordinates (r,@, x), we set 
B = (B,,0, B,), 2/20 = 0, and v = (0,0,v), 


7 9, 
where v is constant. The velocity condition may be used as an approxi- 
mation when there are small perturbations on a uniform flow. With the 
above assumptions E could only have r and x components under special 
circumstances. Considering the line integral of E at constant r and 2, it 
follows thet Z, would be non-zero only if there were a rate of change of 
flux through the area enclosed. We assume that @/ét = 0 everywhere, 
and E = 0 in the fluid. Clearly the applied currents must be such that 
@j,/€0 = 0; we shall also find that j, can only have a 6 component. Ideally 
applied currents should form closed loops, but as in classical electro- 
magnetism, the theory is intended to serve as an approximation for situa- 
tions where the applied current wires take the form of coils. 


Case A(ii). In cartesian coordinates (x,y,z), we set B = (B,, B,,0), 


have v = (v,0,0), where v is constant, but, if the fluid is incompressible, 
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any plane potential flow is permissible. We also make the restriction 
that E = 0 in vhe fluid, and assume that there can be a return path for 
both the applied currents and the currents within the fluid. If case A (ii) 
is regarded as the limit of case A (i) as r + 00, Z, would necessarily be zero. 

Case B(i). In cylindrical coordinates B = (0, By, 0), 2/26 = 0, and the 
velocity conditions are similar to case A (i). 

Case B(ii). In cartesian coordinates B = (0,0, B,), @/éz = 0, and the 
velocity conditions are similar to case A (ii). 

The plane magnetic field conditions for the cases A and B are the same 
as those assumed by Zhigulev (2) when defining boundary layers of the 
first and second kind. 

It is well known that, if j, = 0, the magnetic force term, j x B, in the 
equations of motion is adequately described for case B (ii) in terms of the 
gradient of a scalar magnetic pressure, B?/2u. It follows that vorticity is 
convected along streamlines in incompressible flow when viscosity is 
neglected. Kemp and Petschek (3) have discussed the creation of vorticity 
at a point where the flow crosses a current sheet, i.e. j, # 0. If there is 
no current sheet crossed by the flow, there is no creation of vorticity, and 
we have the possibility of a potential flow solution which is valid for all 
field strengths. The kinematic solution for the incompressible flow condi- 
tions of case B(ii) can then be dynamically exact. It can also be shown 
that vorticity is not created when the normal component of jy is zero at 
a current sheet. However, in potential case B(ii) the fluid pressure 
decreases as B, is increased; there is a greater tendency towards cavitation. 


2. The kinematic equations 
For case A fields it is convenient to use the vector potential A, where 
curlA = B, divA = 0. (6) 
For case A(i), A = (0,A9,0), and for case A(ii), A = (0,0,A,). The 
second equation of (6) is necessarily satisfied. Then under the conditions 
for case A equations (1) to (6) give 
V?Ak = po(v.V)Ak—pijp, (7) 
where A is the magnitude of the vector potential, and k is a unit vector 
in the @ direction in case A(i) and in the z direction for case A(ii). At 
any point j, can only have a component in the k direction. For rectilinear 
flow we introduce non-dimensional coordinates given by 


R=jypourr, X= fyorr, Y = }pory. (8) 
é ( aA Oke _ 9@Ag_ 4Hion 


@R\R @R eX? “eX (pov)? 


(7) becomes 








(9a) 
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A, @A, 404, Siti 


for case A (i), and @x* oY? “ax (pov)? 


- (9b) 
for case A (ii). 

For the plane potential flow of an incompressible fluid we may use 
a method which Boussinesq (4) applied to the problem of heat convection. 
The analogy between the kinematic problems for the case A (ii) and those 
of heat convection in plane potential flow has been shown by Leigh and 
Sutton (5). We define X+i¥ = f(x+iy), (10) 


where grad X = }yov. It is assumed that f is analytic. In terms of the 
coordinates (X,Y), (7) becomes 


@A, @A 902A, . {faX\? | (aX\*)- 

i cee eee SE el i Mars eons ; 1] 

ax? oy! “ax mie (Zz) +(F) j mn 

Thus (9b) may be regarded as the equation for plane incompressible flow 

provided that v is now taken as the local velocity magnitude, v(X, Y). 

Having introduced the transformation to new coordinates, it is convenient 

to return to the use of vector notation, where the operator, V?, is given by 
a oe 

= axat ays 

for case A(ii), and, if Y and Z are cartesian coordinates such that 

Y?+Z* = R*, 


2 


a2 22 a2 

eh hg! CME Lah 

“a ax? ays aga 
for case A(i); (9a) and (9b) can then be written 


(V2—1)e-* Ak = ee. (12) 
(pov)? 
For case B fields the direction of B has already been assumed, and (3) 
is necessarily satisfied, We may find the magnitude of the field strength 
for case B conditions from 


V?Bk = po(v.V)Bk—p curl j,. (13) 
Thus at any point B can have no component in the same direction as k 


for either case B (i) or case B (ii) conditions. Using the coordinate systems 
of (8) and (10), we find 


—2e-X 
(V8 1)e-7 Bk = — th 
pov 
where the operator curl, as well as V?, is applied in the transformed system 
of coordinates, i.e. the magnitude of curl j, is 
or _ Sox 
oX OY 


(14) 


for case B (ii). 





ON SOME KINEMATIC PROBLEMS 323 


The derivation of the case A and case B kinematic equations has not 
been given in full since it differs from standard forms only in the applied 
current density term. A less common equation in magnetohydrodynamics 
is that for the scalar potential of the electric field, U, where 

grad U = —E. (15) 
Taking the divergence of (5), and using the conditions of zero vorticity 
and constant o, we obtain 


V?U + (curl B).v = div jo, 
o 
Using (1), (4), and (5), and eliminating curl B, we have 
V2U = pov.VU+ div Jo— HOV Jo. (16) 
Co 


Since U is scalar, equation (16) is easily applied to field geometries other 
than case B; the only assumptions are those of zero vorticity and constant 
o. However, we shall apply (16) only to case B conditions here. Using 
the coordinate systems of (8) and (10), we obtain 


y (div jo 2dox), (17) 


where the operators V? and div are applied in the transformed coordinate 
system. 


0,-xX 
(V8 1)e-FU = 
a(jov 


3. Boundary conditions 


If the flow field is bounded by a stationary medium, the magnetic field 
should be matched to a solution of the usual magnetic field equations in 
that medium. There are, however, special cases for which the trectment 
is simple. 

With case A fields an iron boundary imposes the condition there of zero 
tangential field component. A superconductor, in which the magnetic 
field is zero, would impose the condition of zero normal field component. 
For the short duration of many ionized gas experiments copper may be 
used to approximate a superconducting boundary condition. 

For a case B field bounded by an insulating medium we have, from 
curl B = 0, that B, is constant in plane geometry, and By is inversely 
proportional to r in axial symmetry. It should be noted that contact 
resistance between metal boundaries and the flow inhibits the passage of 
currents, and metal can be often regarded as fully resistant. If good 
contact between a high conductivity metal and the flow were to be made, 
we would have the condition of zero tangential j component. 

The plane solutions derived below will be for an unbounded flow, but 
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it will be seen that under certain circumstances the classical method of 
images may be used to satisfy the conditions discussed above. 


4. General solution 

The case A and case B kinematic equations are valid only for the plane 
and axially symmetric conditions given above. However, in order to find 
a solution to the equations which can be applied equally to both types 
of geometry, we shall at first consider the following general equation in 
three dimensions: (V2—1)F =f, (18) 


where F and f are unknown and known scalar functions respectively. 
(18) bears the same kind of relation to our type of vector kinematic 
equations as Poisson’s equation bears to the classical equations for mag- 
netic fields in free space. A solution of 
(V?—1)F = 0, 
is e~?/p, where p is the radius in spherical polar coordinates. Let us suppose 
that a region 7 is bounded by a small sphere ©, and a large sphere L,, both 
of which are centred at the origin. The radii of the spheres are p, and p,. 
Green's theorem gives 
f {CPE eM ass f (Cee el as 
p } ép ép | \\ p J ép ep | 
ras [ (e* v2.1) F— F(V?2— y—| a 
1 lp p 
where dS is an element of area on the surface of a sphere. If dQ is the 
solid angle subtended by an area dS, 


—p,e-P | OF 10—e-P(p, +1) | F dQ+p,e-* | OF 0+ 
cp cp 
xX >>) 


+e-A(e+1) { Fly wd of de. 
p 
Xs 


T 


Assuming that @F'/ép is finite at p = 0, @F/ép and F are finite as p > , 
we obtain when p, > 0 and p, > 


1 fe? 
ay Ped f° 9 
F = | ~ fae (19) 


Ps 
T 


Suppose that the case A kinematic equation (12) had no restrictions 
placed on the vectors, and could be applied generally in the cartesian 
system defined by the coordinates (X,Y, Z). Then, using (19), 
pex ( 4e-P-X; . 


A =i | —___j..» - dr, (20) 
X.Y.Z dn | (uav)tp 2°X-F 2 
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where the subscript j is used on coordinates which refer to the position 
of an element carrying current density jy. p will be given by 


p? = (X—X,*+ (¥—¥,*4+(Z—Z,). 


If j, has the appropriate form for a case A (i) magnetic field, (20) represents 
a valid solution to the kinematic problem when X has the definition given 
by (8), and we define Y and Z so that Y*+ Z? = R*. Similarly, (20) is 
also valid for case A (ii) magnetic fields when j, is appropriate, and X and 
Y are defined either by (8) or (10). The fact that Z coordinates have no 
physical meaning in plane potential flow of an incompressible fluid is 
immaterial; we assume ¢/0Z = 0, and the volume integral will extend 
from Z; = —« to Z; = . A single line current I, is given when j, acts 
on an infinitesimal area 5S in the (X,Y, Z) system, and 


4, 58 
(uov)? 85-0 


x —p-X; 
(20) reduces to Ant. I, dl, 
4 p 
where di is measured on the line current in the (X, Y, Z) system. 
Similarly, the solution of (14) is 


X f 2e-p-X 
S a {= — eurl j, dr. 
4n J poup 


It is convenient to express the solution in terms of sheet currents. We 
shall restrict the analysis to case B (ii) although similar methods could be 
used for case B(i). Suppose that the coordinates m and n refer to direc- 
tions parallel and perpendicular to j, respectively, where (m,n,z) is a 
right-handed system. For a sudden rise in the magnitude of j, from zero 
to Aj, we have curl jo| 8n =—> Ajo. 
A sheet current with magnitude J, is given by 
2 . 5 J 
= Ajo dn aap Io 


The solution for Bz due to a small element carrying a current density 
Aj, is 


x5 p9 —p-X 
B, = pe a ile dmdZ ; n. 


pavp 


4m On | 
—@ 


Reversing the order of differentiation and integration, and noting that 
the physical length of the element 28m/yov should be constant with 
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respect to the differentiation, we obtain 


_ mek f 0 fer% 
B, = = fat - )io 8maZz, 


The solution for the electric scalar potential in equation (17) can be 
taken directly from (19), giving 


= ex 2e-P-x 
U = pee *e- ’ (divj,— 2jox) dr. 
By a similar method to that used for magnetic field of case B solution, we 
obtain for the plane case 
x Qe-P- h é fe-P 
Dm of a and, + — f 2h, 2: 
i: | 5 ox bm += = | =| : —h SmdZ,. (23) 


—~2 —@ 


5. Some particular case A solutions 


The simplest case A field is that given by a straight line current. 
Taking the line current to be along the Z-axis, equation (21) gives 


Az — Ho .x [ 


e 4X*+¥*+Zpt 
an) (XE+¥EC ZF 
Substituting + Z, = (X*+ Y*)'sinhy, we obtain the standard form 


~ 
ul, .. a6. ati 
A cin, ge et X*+¥ *ooshd dy, 
Zz Qa p 


0 
which gives (see for example McLachlan (6)) 


Ay = K%eXKy(X*+Y2), (24) 


where A, is a modified Bessel function of the second kind and of zero 
order. (24) is similar to the solution found by King (7) for the temperature 
due to a heat source in potential flow. The magnetohydrodynamic result 
has been derived by Falk (8) and by Ludloff (9) for rectilinear flows only. 
For small values of X and Y 


y oe lon x? + y2)t, 


The magnetic field near the wire is approximately the same as in free 
space, and since the field is hardly distorted, large velocity perturbations 
close to a slender wire will have little effect on the solution at moderate 
X and Y. 
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The potential flow past an iron wall with a semicircular projection is 
given by the well-known Joukowski transformation, 


2 

X+iY = }yor, (r+iy+ | (25) 
where v, is the free stream velocity, and d is the diameter of the semicircle. 
The streamline which follows the boundary of the wall is given by Y = 0. 
If we suppose a line current at (X,,¥;), the boundary condition at the 
iron can be satisfied by taking an equal image current at (X,, —Y;), and 
the field is 


Az = MeX-XAK(X—X))2+ (YY P+ KX XPV +H). (26) 
This example shows the limitations of the method of solution. Even for 
simple boundary conditions we are restricted to problems where an image 
method is possible in the (X,Y) plane. Equation (26), where X and Y are 
given by (25), is also the solution for symmetrically placed line currents 
in the flow about an iron cylinder, but the solution for an arbitrary 
distribution of line currents is not easily found. The particular case of 
line currents which lie at the stagnation points of a cylinder has been 
calculated for values of the magnetic Reynolds’ number, Ry = povgd, 
equal to 2 and 10. The field lines are shown in Figs. 1 (a) and 1(6). The 
line currents are of equal magnitude, but opposite in direction, and the 
number ascribed to each field line is the value there of 47A/uJ,. For the 
higher magnetic Reynolds number the field lines are swept strongly down- 
stredm; the cylinder begins to form a magnetic wake. For very high 
magnetic Reynolds numbers the field is confined to a thin boundary layer 
on the cylinder and to the wake behind the cylinder. 

An important solution for axial symmetry is that for the field due to 
a wire loop immersed in a uniform stream. If in the non-dimensional 
frame the position of an element of the loop is given by (Ro, 6;,0), where 
R, is constant, (21) gives 





= et Re+R*+ X*-2RR,cos 6, : 
Ag = Hg j Cr Rex? 7 SRR, condi 0" 6; d0;. (27) 


re 


If R®4+X* > R?, and if the magnetic Reynolds number of the loop is 
small, so that Ry < 1, (27) gives approximately 


plo Ry X~+R*+ xt ge 
Ae = Se (REE XIN es 


RR, cos 6; if 4 1 
( R2+- XI (R?4-X?)t 





}}eos 0, d8,, 
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Free stream -1-0 
— 


Fie. 1 (a). Field lines about an iron cylinder with line currents at the 
stagnation points—Ry = 2. 


Free stream —0-25 
——— 


Fic. 1 (6). Field lines about an iron cylinder—R,y = 10. 





2 
that is, Ag = Hig Se eX 4 : R (28) 


(R24 X?)t) R24 X?" 
The solution (28) is that for a magnetic dipole at the origin. 

The magnetic field on the axis of a finite loop is easily evaluated. The 
result is of no great significance in magnetohydrodynamics if we are using 
the kinematic problem as a first-order approximation for finding magnetic 
forces; by symmetry there can be no force on the flow along the axis of 
the loop. Here we shall use the field strength on the axis of the loop as 


2r 
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a guide to how much the magnetic field is distorted near a coil. As R > 0, 
it can be shown that 


(29) 





Bx l Ro 

“ply (B+ X28] RE+X* 
The value of By r,/zJ, is plotted against X/R, in Fig. 2 for values of Ry,, 
based on the diameter of the loop, equal to zero, 4, and 20. As in the 
example of the iron cylinder, the field perturbation can be explained by 
imagining the field lines swept downstream. The point of maximum field 


= fet mexH(1 4 


Byro 
al 


0S, 





Fic. 2. Field on the axis of a wire loop. 


strength occurs at X/R, > 0, if Ry > 0, but, as Ry, is increased, the 
penetration of the field from the loop into the flow is reduced. 


6. Some particular case B solutions 

For the plane case B conditions we shall treat only three basic problems 
from which it will be seen that more general solutions can be constructed. 
The field due to a small element of a sheet current path can be found from 


22) as 


B, = — sex [ (74a p)—l/phe-P- BY 42+ 
“i p 


4 Ho x T—Ha4y p)je?-*s 8X, dZ,, 
4n p* ; 


-~@ 
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where 5Y; is the extent of the element in the Y-direction, and 5X ; is the 
extent in the X-direction. Assuming that the element lies on the Z-axis, 
and making the substitution + Z, = (X*+Y*)tsinhy, we obtain 


= MeN xk (X24 ¥8)— 


(1+ ies 
rryveke'+ acres 








Ae) 


le —(X*+¥F% coshy dis uh 








+e x Y (14 1 le 4(X*+¥Ft cosh yd db. 
J (X#+Y?)t cosh y\ (X2+ ¥2)* cosh & 


Integrating by sr terms of the type 
e+ X*+F Yt coshy 
cosh%, =” 
the remaining integrals reduce to standard forms, and 
eee x 
By = MoU teX{K(X*+¥4)— aaaall yy + 





+ Kyixr+¥7y, (30) 


“lara 
An isolated element of the sheet current is in fact an electrode a 
and (30) may be taken in two parts, the magnetic field due to an electrode 
doublet whose axis is perpendicular to the flow, and the field due to a 
doublet whose axis is parallel to the flow. Current streamlines, for 
values of 47Bz,/yJ, 5m = constant, are shown in Figs. 3 and 4. A case A 
magnetic doublet would give a solution for Az in the same form as (30). 
Figs. 3 and 4 could thus be interpreted also as the field lines due to 
magnetic doublets. 

The scalar potential for electrode doublets is found from (23), giving 
when the axis is perpendicular to the flow 


J 8Y; Y rs\a\ 
ie sotet espa Kult +. (31) 
The equipotentials have the same form as the current streamlines for 
a doublet whose axis is parallel to the flow. The second doublet solution 
gives 

ee J, 8X; bX, | xX 
“2ra- (X22) 

We can form a final basic solution for case B(ii) from the doublet 
solutions. Suppose that there is a doublet at (0, Y;) with J, lying in the 
X direction, and that there is a second doublet at (0, Y;+-8Y;) with J, in 


U 


2ro 


K(X*+¥%)+ —~__K,(X*+¥2'}, (32) 
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xr=0 05 1-0 
t 1 — 





Fic. 3. Current streamlines for an electrode doublet, axis 
perpendicular to the flow. 


x=0 0-5 1-0 
— | 


Fic. 4. Current streamlines for an electrode doublet, axis 
parallel to the flow. 





332 M. D. COWLEY 


the —X direction, i.e. there is an anticlockwise circulation of J,. From 
(30) we obtain 





bX, @ (Y—Y,) ae 
BD _ wh Sees j 2 —Y,)2)t f 
oo a leree rp + O—2 a, 
After differentiating, Y, is put equal to zero to give the solution for two 
opposed doublets at the origin, 


5X, dy, { 1 ‘ y2 
B, — te j°*i eX 24 ys) __"____se FKy( X24 y2)t}, (35 
Zz Qn e xe yap Mal X +Y ) (xe yas MelX +Y ) (33) 
where the recurrence formulae for modified Bessel functions of the second 
kind have been used. Similarly for two opposed doublets whose axes are 


perpendicular to the flow 





2X+1 ,, , 
(xa ya AXP +4 + 
2 


B, ai tot ex — Ky X?+ Y2)t+ 
7 


x vin oe 
— apy erry, (34) 


where the directions of J, have been again chosen to give an anti-clockwise 
circulation. Adding the two results, we obtain 


p, — we®X;8Y; x 
7 





oF : pay Kl XP+ ¥9))— Ky X84 ry). (35) 
Since the current sheets have been arranged to form a closed loop, (35) 
is the solution for the magnetic field due to an infinitesimal solenoid. 
A region filled by infinitesimal solenoids, each carrying density J), is 
equivalent to a finite solenoid, current density J,, bounding the region. 
The form of the current streamlines given by (35) is similar to that for 
an electrode doublet whose axis is perpendicular to the flow. In finding 
the magnetic field for a finite solenoid, the constant Jj, must be added 
to the solution for B, within the solenoid. 

Figs. 3 and 4 show that the current streamlines are swept downstream 
by the flow. For the doublets a particular value of B, is now found to 
occur downstream of the point where it would be found in a stationary 
medium. In contrast the case A fields are such that lines of constant A 
are swept downstream, and the variation of B depends on the geometry. 
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SUMMARY 

Exact solutions of the equations of magnetohydrodynamics are shown to exist 
for the radial flow of a viscous incompressible fluid between non-parallel, plane 
walls. Both diverging and converging flows are discussed, using the approximation 
of small magnetic Reynolds number. It is found that the maximum permissible 
channel angle for purely divergent flow can be increased without limit by a suffi- 
ciently strong azimuthal magnetic field. Solutions of boundary layer type are 
obtained for converging flow at high Reynolds numbers, and it is shown that the 
magnetic field reduces the displacement thickness of the boundary layer. 


1. Introduction 

THE equations of magnetohydrodynamics can be solved exactly for the 
case of two-dimensional steady flow between non-parallel plane walls, of 
a viscous, incompressible, electrically-conducting fluid; this is a straight- 


forward extension of the celebrated Jeffrey-Hamel problem in ordinary 
hydrodynamics. In the present paper we examine approximations to some 
of these solutions, which are valid for small values of the magnetic Reynolds 
number. The problem is of considerable theoretical interest, particularly 
as solutions can be obtained in closed form; there may also be practical 
applications of our results to magnetohydrodynamic flows in nozzles and 
diffusers. 

A survey of the Jeffrey-Hamel problem has been given by Goldstein 
(1, 105); in this paper we follow Goldstein’s treatment closely, and the 
reader need not refer to any of the earlier literature. However, mention 
must be made of the work of Rosenhead (2), who treated the problem 
more fully than earlier writers, and Millsaps and Pohlhausen (3), who 
considered thermal effects; there are also surveys in the books by Landau 
and Lifshitz (4), and Dryden, Murnaghan, and Bateman (5), which are 
rather more extensive than that of Goldstein. 

In section 2 of this paper we show that for the Jeffrey-Hamel problem, 
the equations of magnetohydrodynamics can be reduced to a set of three 
ordinary differential equations, two of which are linear and of first order. 

+ Now at Defence Research Telecommunications Establishment, Shirley’s Bay, Ottawa. 


{Quart. Journ. Mech. and Applied Math., Vol. XIV, Pt. 3, 1961] 
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It is necessary that the magnetic field should originate in a combination 
of a line magnetic pole and a current, both coincident with the line of 
intersection of the channel walls (the z-axis in cylindrical polar coordinates). 
The boundary conditions to be satisfied at the walls are also given in this 
section. 

One of the features of bounded flows in magnetohydrodynamics, is that 
the flow causes changes of the magnetic field in the solid regions, and as 
a result the solution for the field must cover all space if the problem is 
to be treated properly. In section 3 we examine this aspect of the problem, 
and find that for all conditions to be satisfied, the total flux of fluid from 
the z-axis must be zero. Thus purely converging or diverging flow can 
only be obtained in a channel if there are other channels which allow the 
net flux to vanish. It is found that the boundary conditions in the exact 
problem (with a given configuration of channels and solid regions) are not 
at all simple, due to the interaction of the flow with the external magnetic 
field. With the approximation of small magnetic Reyno!ds number how- 
ever, the boundary conditions and the equations take on their simplest 
form, provided also that there is no magnetic pole on the z-axis. This 
result, which is obtained in section 4, is most fortunate since the magnetic 
Reynolds numbers obtainable in practice are invariably small. 

We consider purely diverging flow in section 5, following the treatment 
given by Goldstein for the non-magnetic case. It is found that the mag- 
netic field opposes the formation of regions of back-flow, which would 
otherwise always appear at sufficiently large Reynolds numbers whatever 
the angle of divergence of the channel; in fact back-flow can be completely 
prevented provided 3R,,S > 2, where R, is the magnetic Reynolds 
number and S the Alfvén number. In the case of converging flow at large 
Reynolds numbers, which is examined in section 6, the magnetic field 
reduces the boundary layer displacement thickness and increases the 
curvature of the velocity profile in the boundary layer. Thus we find 
Cowling’s (6) description of the effect of a transverse field, as a sort of 
‘viscosity’ tending to destroy motion across the lines of force, to be slightly 
misleading; in fact the magnetic field operates in exactly the opposite 
manner to viscosity since it tends to prevent the diffusion of vorticity 
and the appearance of back-flow. 


2. Equations and boundary conditions 

Consider a system of cylindrical polar coordinates (r,@,z) in which the 
channel walls lie in planes 6 = constant, and intersect in the axis of z. 
It is assumed that there are no changes with respect to z, that the motion 
is purely radial and that there is no magnetic field in the z-direction. 
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© (ru) = < (rH) +S 
cr 


fF. 43.  A6 
where V2 == +2 xa) 
v is the kinematic viscosity of the fluid, » the magnetic diffusivity, p the 
density, and » the magnetic permeability. The velocity and magnetic 
field are of the form (u, 0,0) and (H,, Hg,0) respectively, and y represents 
the ‘total’ pressure; thus if p is the fluid pressure, then 


x = p+ £(H}+H}). (2.7) 


The equation of continuity (2.5) requires that u should have the form 

u = Uf(0)/r, (2.8) 

where U is a positive constant having the dimensions (velocity x length). 
To match this simple dependence on r, we must choose the magnetic field 
such that H, = Hg(0)ir, Hp = Hr, (2.9) 


where H is another constant which is determined by the strength of an 
isolated current flowing along the z-axis. It is not necessary that H, should 
be zero when Hy == 0, but as the magnetic field has no effect when it is 
purely radial we shall ignore this possibility and use the same H for both 
components of the magnetic field. The form for x consistent with (2.8) 
and (2.9) is found to be 

x = ph(@)/r?+constant. (2.10) 
Substituting these forms in equations (2.1) and (2.2), we obtain two 
ordinary differential equations, as follows: 


f2?+2h/U +of"/U +(g' —g?—1)uH?2/(4apU2) = 0, (2.11) 
h’ = 2vUf"’, (2.12) 
dashes being used to denote differentiations with respect to 6. Equation 


+ See, for example, Yih (7). 
Zz 
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(2.12) can be integrated at once, giving 
h = WUS+C,, (2.13) 
where we use C, to denote arbitrary constants as they occur. Reynolds 
and Alfvén numbers are defined ast 
R= Uj, S = pH?/4npU?, (2.14) 
respectively, so that (2.11) becomes 
f?+4f/R+f"/R+ S(g'—g?—1)+C, = 9, (2.15) 
on substituting (2.13) and putting C, = 2C,/U?. 
As is commonly found in magnetohydrodynamic problems of this type, 


equation (2.3) reduces to the derivative of (2.4), the latter equation 


becoming g' +Ryf = 0, (2.16) 


where we have introduced the magnetic Reynolds numbert 
Ry = U/». (2.17) 
Equation (2.16) is simply Ohm’s law, and it implies that there is no electric 
field in the z-direction; hence the total electric field is zero. The axial line 
current which produces the azimuthal magnetic field is assumed to flow 
in an effectively perfect conductor, thus no electric field is required to 
maintain it. 
The boundary conditions appropriate to this problem are as follows 
(taking the channel walls to be 6 = a,): 
(a) There is no slip at the channel walls, 
fla) = 0. (2.18) 
If only symmetric solutions are considered, 
f'(0) = 9, (2.19) 
6 = 0, being the midpoint of the channel; whatever the case, df/d@ is zero 
at some point in the flow, since f = 0 at the walls. 
(b) The tangential component of the magnetic field must be continuous at 


the channel walls; thus if H, = H’g,/r within the solid boundary at @ = a, 


say, then Hg(a,) = 9, H’. (2.20) 


Applying this condition at the other wall (@ = a,) we see from (2.16) that 
H, = H’g,/r within the solid boundary at 6 = a, where 


Xe 
H'(g,—92)/Ry H = | f d@ = total radial flux of fluid in the channel. 

“ (2.21) 
+ For the purposes of this paper these definitions of R and Ry are quite satisfactory. 


In general, however, Q should replace [’, where Q is the volume flow of fluid in the channel 
fe.g. (2, 4)}. 
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Hence, as far as the solid material on either side of the channel is con- 
cerned, the region of flow appears as a current sheet which provides a 
jump in the tangential component of the magnetic field. 


(c) The normal component of the magnetic induction at the channel walls 
must also be continuous. Thus if Hy = H'/r in the solid, then 


pH wt’, 
where yu’ is the magnetic permeability of the solid. 


(d) The conditions involving the electric field are satisfied automatically, 
as it is zero € verywhere ; 


It is noteworthy that the conductivity of the solid material does not 
affect the problem in any way. 

The above conditions will not be sufficient to provide a unique solution 
for every set of values of the parameters involved; in fact it is found in 
the ordinary hydrodynamical problem that eigensolutions are obtained 
for the velocity distribution. Presumably this behaviour will occur in 
the magnetohydrodynamic problem also, but as we are interested only in 
certain special aspects of the problem in this paper, we need not be con- 
cerned with a possible lack of uniqueness in general. 


3. The solution for the solid regions 


When magnetohydrodynamic effects are included in the Jeffrey-Hamel 
problem, it is not sufficient to consider the motion of the fluid as being 
due to a line source or sink coincident with the z-axis; in fact it is shown 
below that the net flux of fluid from the axis must be zero. With this 
restriction it is clear that we cannot examine all types of flow unless there 
is more than one channel; such a situation is shown in Fig. 1. In the 
diagram and in the working which follows only two channels are con- 
sidered with walls at @ = a,, a, and @ = ay, a, respectively. However, 
the generalization te allow for any number of channels is trivial. 

In regions B and D of Fig. | let H, have the values H’g,/r and H’gp/r, 
respectively; then from condition (b) of the previous section we have 


Hg(a,) = Hg(as) = H'gp (3.1) 
and Hg(a,) = Hg(a) = H’g,. (3.2) 
Similarly, from equation (2.21) we see that 

g(a)—gla,) = Ry | f de (3.3) 


™ 


g(%4)—glas) = Ry | f a8. (3.4) 


a3 
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the constants U and H being the same for each channel. Combining the 
above results we obtain 


ae 


J 440+ [sao =o. (3.5) 


which shows that the total radial flux of fluid must be zero. 


2 


Fic. 1. Diagram of the (r, @) plane for a configuration with two channels 
(A and C) bounded by walls at @ = a,, a, and 0 = as, a, respectively. The 
hatching denotes the boundaries of the solid regions B and D. The angles 
subtended by the four regions A, B, C, and D are 27a, 2mb, 2nc, and 2nd 

respectively, where a+6+c+d = 1. 


The current-density induced by the motion of the fluid across the 
magnetic field is in the axial direction and is given by 


iflé 1 oH, ee a Ry H 
=|; 5¢%)- : |= — 7590) = *Sf(0). (3.6) 


r or r 0 


The current flowing in channel A between r = r, and r = r is therefore 


=u= tf (2 aaar - = af ff na (3.7) 


Tr, Gy 
and the total current flowing in both channels is 


im “4, Ay j f(0) d0-+. f f(6) a} = 0 


r3i—70 4n 
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whatever the value of r. Hence there is no accumulation of charge at 
infinity due to unclosed currents; currents induced by outward flowing 
fluid are exactly balanced by currents in the opposite direction induced 
by inward flowing fluid. This leads to a further result which relates the 
azimuthal field to the strength, 7, of a current flowing in a wire coincident 
with the z-axis; for 


[ H,r dé = total current flowing within the contour. (3.9) 


4m 


r 
Thus = }H(a+c)+H'(b+d)}. (3.10) 
Applying Gauss’s theorem to the cylinder generated by the contour I, 
we obtain $ . 
| pH,r dd = | wHg(6) dé = 2nM, (3.11) 
r r 
where M is the strength of a line magnetic pole lying along the z-axis. 
Using (3.1—4) and (2.22), equation (3.11) can be written 
27M . “ « ‘ 
wr 2a(bgx+dgp)+ | g(6) d0+ | g(9) dé; (3.12) 


integrating by parts and substituting for g’(@) from (2.16), this becomes 


a a 
27M * a, 


_—- 2n(a+b)gp+2n(c+d)gpn+Ry | Of(0) d+ Ry | Of (0) dé. 
_ i (3.13) 
The integrals in this last equation represent changes in the distribution of 
the radial field due to the motion of the fluid. If there was no motion 
these would vanish and the magnetic field would be given by 


In = Ip = M/pH. 


It is noteworthy that the changes produced by the motion are O( R,,) 
and in particular, if M = 0, 


2n(a+6)gp+2n(c+d)gp = —Ry| | Of (8) d+ | 6f(9) i}. (3.14) 
Combining (3.14) with (3.1—4), we see that both g, and g, are O( Ry) 
when M = 0. This result is of some importance when the magnetic 
Reynolds number is very small, for only then do the boundary conditions 
for g(@) take on a manageable form. 


4. The approximation of small magnetic Reynolds numbers 
Although one can, in principle, solve the magnetohydrodynamic Jeffrey— 

Hamel problem exactly, it is of interest to examine the simpler cases and 

in particular those which have greatest relevance to practical conditions. 
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Thus, bearing in mind the small magnetic Reynolds numbers attainable 
in practice, we will consider an approximation which is valid in these 
circumstances. As it happens, this is also the case for which the equations 
and boundary conditions take on their simplest forms. In general, despite 
the simple relationship between g’(#) and f(6), the problem requires the 
solution of a non-linear, third-order differential equation (found by putting 
f(@) = F'(6) and g(6) = —Ry F(6) in (2.15)) with some very awkward 
boundary conditions. However, if R,, <1, we have 


g(9) = 9, +O( Ry), (4.1) 
and so l—g’+g* = 1+ Ry f+g?+O(g, Ry). (4.2) 
Hence, if g, = O( Ry), we can ignore the term O( R4,), so that the equation 


to be solved is 
r+(s os S\f+3f" Ne (4.3) 


where we have written C; = C,—S. It will be noted that this equation 
is of second order instead of third order, and that the reduction has been 
effected without neglecting high order derivatives. Equation (4.3) is 
similar to the equation found in the ordinary Jeffrey-Hamel problem, the 
only difference being in the modification of two of the coefficients. We 
might have obtained the equation by an approach similar to that of 
Rossow (8), but the present method is more rigorous and emphasizes the 
necessity for having M = 0. 

Multiplying equation (4.3) by 2f'(@), we see that it can be integrated 


to give 1 /df\? 
alas) = C,—20,f— (p—Ra s)f* —T, (4.4) 
which can be written in the form 
2 2R 
(ih) =e MDD) = AST), say. (4.5) 


The roots of the cubic T(f) (viz. f,, f., fs) must satisfy the following 
relat ear : 


hhh =! lethal =n frthet hs = 1h 8-5 


(4.6a, b,c) 
fi, fe, and f, may all be real, in which case we choose f, > f, > f,; other- 
wise one root (f,) is real and the remaining two are complex conjugates. 
Since f’? > 0, all values of f occurring in the flow lie in a connected region 
in the T(f) diagram with T(f) > 0. The no-slip condition requires that 
f = 0 is in this region; and as f’ is zero at one or more points in the flow, 
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then f = f,, fz, or f, at such points and the region must also include at 
least one of these values of f. 


5. Purely divergent flow 

In the non-magnetic problem, it is found that purely divergent flow is 
possible only when the angle between the channel walls is less than a 
certain value which decreases with increasing Reynolds number. The effect 
of a magnetic field on the relationship between this maximum permissible 
channel angle and the Reynolds number is examined in this section. 

For purely divergent flow f(@) is everywhere positive or zero and thus 
the solution is defined by a connected region in the 7(f) diagram with 
T(f) > 0 and f ranging between zero and a positive root of T(f) = 0. 
There are three cases in which this is possible; these are called types (a), 
(b), and (c), and are shown schematically in Figs. 2, 3, and 4. It is only 
necessary to consider symmetric solutions and so it is convenient to take 
the walls to be at @ = +a. Hence the upper limit (f,, say) of the range 
of values of f corresponds to the flow at the midpoint of the channel 
(@ = 0); similarly condition (a) shows that the lower limit (f = 0) corre- 
sponds to the walls. Between these limits, 7'(f) > 0 for types (a) and (6), 
and so df/d@ cannot change sign; thus 

df 2R 
B= = [7 \runy. (5.1) 


with the upper sign referring to 0 > @ > —a, and the lower sign to 
x > 06> 0. The solutions can therefore be written in the form 


a—|0| = M&) fren tf, (5.2) 


fm 
: 3\ fF , 
which gives aq = (ez) | [7(f)|-* df. (5.3) 
0 


Slight modifications of these results are necessary for type (c) flows. 


Type (a): Here f = f, > 0 at the midpoint of the channel, and it is 
convenient to choose U so that f, = 1. Since f, and f, are negative, zero, 
or complex conjugates, then f,f, > 0 and we see from (4.6a) that C, > 0 
also. Substituting for (f,+f,) from (4.6), 


Tf) = (\—f)|8C— (BR S—F— 1) +9". (5.4) 


In order to make the integral in (5.3) large, we must ensure that 7'(f) is 
as small as possible in the range concerned. This is done by taking C, = 0 
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so that 7'(f) has a zero at each end of the range, and the integrand becomes 


infinite at these points. Thus « is a maximum when f,=0 and 
f, = $Ry S—6/R—1, as shown in Fig. 2.2. Provided 3R, S < 6/R+1, 


1 

the integral f [7(f)|-* df is convergent, as the singularities of the inte- 
0 

grand are due to factors f- and (1—f)-+. Thus we obtain a finite value 


Of Omax 1 
amex = ,/ (52) | [an] s+ (143 —Bears)z]) ay (5.5) 


TY) T(f) TS) 


Alam, WE NIB», 


1 2 3 
Fic. 2. T(f) diagrams for purely divergent flow, type (a). The hatched regions define 
the values of f which occur in the flow. (2.1) Typical case; /, and f, are negative, 
zero, or complex conjugates. (2.2) Optimum case for given R and Ry, S, resulting in 
the maximum possible value of « for purely divergent flow; f, = 0, f, < 0. (2.3) The 
limiting case of type (a); f, = f, = 0. 











If we put f = cos*y, (5.6) 
equation (5.5) reduces to 


toa 
Me eo a iS 
|G) smes <= k | (1—sin*)! = kK(k), (i ./) 
0 


where iP = 2+ 5-1 s] a (5.8) 


and K(x) is the complete elliptic function of the first kind with modulus k. 
The presence of the magnetic field increases the value of k (which would 
otherwise never be more than 0-707) and hence increases the maximum 


permissible value of «. In the limit as k* > 1, the right-hand side of (5.7) 
in 


becomes i secy dy and R'a,,,, increases indefinitely. This limiting case 
0 


is shown in Fig. (2.3). The condition k* = 1 corresponds to f, = 0 so 
that 7(f) has a double zero at f = 0. Thus the integrand of (5.5) is 
(1—f)-*f- and the integral diverges logarithmically. Type (a) flows are 
not possible when $2, S > 6/R-+-1, and in fact, the limiting case, k? = 1, 
is the borderline between type (a) and type (c) flows. 
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Type (6): In this case f = f, > 0 at the midpoint of the channel, and 
it is now convenient to choose U so that f, -= 1. Since f, and f, are both 
positive or complex conjugates, then f,f, > 0. Thus C, > 0 and once 
again 7'(f) is given by (5.4). The maximum value of « is now obtained 
when 7(f) has a double zero at f = 1, as shown in Fig. 3.2, and in fact, 


1 
the integral { [7(f)|-* df diverges logarithmically as in the limiting case 
0 


TUS) T(S) TY) 


\ \ 


2 3 

Fic. 3. T(f) diagrams for purely divergent flow, type (6). (3.1) Typical case; /, 

and /, are positive or complex conjugates. (3.2) Optimum case; f, = f, < f,. (3.3) 
Limiting case of type (6); f; = f. = fs. 














of type (a). Thus purely divergent flow is possible for all values of « if 
fe = 1 and f, = §Ry S—6/R—2>1. The limiting case of type (6) is 
that for which f, = 1 and 7(f) has a triple zero. This is the borderline 


1 
between types (6b) and (c), and it is the only case in which | [T(f)}° df 
0 


diverges more rapidly than logarithmically. 


Type (c): For types (a) and (6), (3R,y, S—6/R—2) is less than (—1), and 
greater than (+1), respectively; type (c) represents the intermediate flows 
with |3R, S—6/R—2| < 1. As in type (5) there is no limit to the values 
of Rta permitting purely divergent flow, for it is possible to have the 
situation sketched in Fig. 4 where 7'(f) has a double zero within the range 


of values of f occurring in the flow. If we take f, = 1 then the position 
of this zero is given by 


f=h=h=4[ws—5-1]. (5.9) 


Since | 7'(f)}* changes sign at f = f,, it is necessary that the alternate sign 
in (5.1) be taken, in order that the maximum permissible value of « should 
be attained; this also ensures that {(@) is a non-increasing function of {6}. 
The appropriate modification to (5.3) is therefore 

, 1 


== /(55) | [T(f)|>\ af. (5.10) 


0 





346 W. I. AXFORD 


It is interesting to note that 7'(f) has zeros at the midpoint of the channel 
and at the wall in type (a) flows, the latter zero being double in the limit- 
ing case. In the intermediate flows (type (c)), this double zero moves away 
from the wall and eventually coincides with 
the zero at mid-channel in the transition to 
type (6). Thus for types (6) and (c), there is 
\ no imminent separation at the walls as there 

is in type (a), and as a result, no upper limit 
on Ria. 

It is clear from the above results that the 
effect of a magnetic field on purely diver- 
$os..6; PLS) dinaneni for purely gent flow is to allow greater channel angles at 
divergent flow, type (c); optimum @ given Reynolds number, than are possible 
case, fs = f, < f;. The limiting without the field; moreover, if Rj, S > 3, 
cases of this type are as shown purely divergent flow is possible at all finite 

ae 5A eae Se. Reynolds numbers without restriction o 
y n 
the magnitude of the channel angle. The variation of R'a,,,, with R for 
given values of R,, 8 is shown in Fig. 5. These curves should be com- 


T(f) 























i i i A iL > FR 
107? 10! 1 10 10? 10? 10° 
Fic. 5. Rta,,, a8 a function of R for various values of Ry S. When 
Ry S < %, the curves become Ric, ~ constant, whereas for Ry S > §, 
the curves become R ~ constant; the former case corresponds to type 
(a) flows, and the latter to types (6) and (c). 








pared with that for the non-magnetic case (Ry S = 0). The asymptotic 
values of Rta,,., for large values of R are plotted against R,, 8 in Fig. 6. 
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Values of R,, S greater than 3 are by no means unreasonable in practical 
applications. Thus our results are of particular interest in magnetohydro- 
dynamic channel flow as they indicate that a sufficiently strong magnetic 
field may prevent separation in diffusers, thereby reducing head losses 
and increasing the efficiency of the channel. For boundary layers in 


{ 


RyS 
0-667 


0-6 











. nh ee ee 
Fic. 6. The asymptotic values of Ria,,,, for large values of R plotted 
against Ry S. 
general, we infer that any tendency to separate is opposed by a transverse 
magnetic field, and may even be completely suppressed by a sufficiently 
strong field. 


6. Convergent flow at high Reynolds numbers 
For convergent flow f < 0 everywhere, and as shown in Fig. 7, T(f) 
must have one root positive or zero, and the remaining two roots real 
and negative; thus f, > 0 and f, < f, < 0. It is convenient to write 
f=, f, = —Ws, f=—vw, (6.1) 
where w,, w,, and w are all positive or zero, and we choose U so that 
te 1. Hence from (4.6 c) 


v 6 J 
wW,—W, = Ry 8— 5+. (6.2) 


The solution can be written in the form 


f (FF) 6\) = { [(w,+-w)(1—w)(w,—w)}-t dw. (6.3) 
0 


There is no restriction on « in this case, but for large Reynolds numbers 
the integral on the right of (6.3) must also be large. This will be true when 
w, > | (so that 7(f) has a double zero at f = —1), and 


=. : ‘ 
w, > §Ry S—7+2 = PRy S+2. (6.4) 
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The solution can now be found easily by substituting 
w = (1+w,)tanh*é—w, 


in the differential equation 


d\6) 
“0 ale eons 


We find that "Ez = -/ | ar 


Ri +w,)} 








la \ 4 


a b 
Fie. 7. T(f) diagrams for converging flow. (7.1) Typical case; all roots 
real, f, > 9, f; < fz < 0. (7.2) Flow at high Reynolds numbers; /, = f, < 0, 
i, >%. 


x 








Hence, using the boundary condition (2.18) in the form 


tanh’ = w,/(l+w,) when (|@| = a, (6.8) 
we obtain 
tanh?é = a i = tanh] {ACE N(a— 0 +7}, (6.9) 
l+w, 6 
w 


where tanh*y = yogis = 1—}(14+4R,, 8)". (6.10) 
1+w, 


Thus 
w = 3(14+-4R,, S)tanh4{{}R(1+4Ry S)}(a—|0|)+y]—2—3Ry 8, 
(6.11) 
which reduces to Pohlhausen’s (1, 143) solution for boundary layer flow 


in a converging channel when R,,S = 0. The dimensionless boundary 
layer displacement thickness, 5, can be obtained directly from (6.6), as 


= 1 
RS = | (1—w) d{ Ra—|0\)} = 3 [ (w,+w)- dw 
0 0 


= 3V2(144R,, S)*—2v3(1+3Ry, 8). 
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Fic. 8. Boundary layer velocity profiles for purely converging flow. The 

lowest curve represents the Pohlhausen profile (RyS = 0), while the 

curves above correspond to Ry S = 0-25, 0-5, 1-0, 2-5, 5-0, 10, and 25, 
respectively. 
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Fic. 9. The variation of the boundary layer displacement thickness with 
Ry S for purely converging flow. When Ry S = 0, Rt3 = 0-778, while 
R16 ~ (Ry S)-* for Ry S > 1. 


The variation of the magnetic field can be found quite easily since 
g = —Ryf = Ryw. (6.13) 
We find after some algebra (using condition (6) of section 2) 
g(9)—g, = Ry G(\0|) when —a <6< 9, 
Jo—G(9) = Ry G(\0|) when a>é>0, 
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where 
G(10}) = a— a [a fatwepttann{ [FEE a) +r} —wt]. 


(6.14) 


For large values of Ri(a— |@|) we see that 


G(\0|) > a— |6|—8, (6.15) 


as one would expect from the definition of the displacement thickness. 
Some velocity profiles for a range of values of R,, S are shown in Fig. 8, 
where it can be seen that the main effect of the magnetic field is to produce 
greater curvature of the profile, and as shown in Fig. 9, to reduce the 
boundary layer displacement thickness. These effects are precisely the 
same as those found by Hartmann (9) for one-dimensional flow between 
parallel plane walls with a transverse magnetic field, and we infer that 
they occur generally in flows of this type. It is likely that the stability 
characteristics of a boundary layer are affected by a transverse magnetic 
field, due to the increase in curvature of the velocity profile. However, 
Lock’s (10) results suggest that the presence of the field should be favour- 
able for the maintenance of laminar flow, although the small disturbance 
theory tends to overestimate the stabilizing effect when R,, S is large. 


ADDENDUM 

Jungelaust has recently discussed some boundary layer problems which are 
essentially solutions of boundary layer type associated with the Jeffrey-Hamel 
problem. 

Since Jungclaus’s symbol N is equivalent to Ry S, it is easily ascertained that the 
solutions for converging flow at high Reynolds numbers are the same (equations 
(17) and (6.11) respectively). The displacement thicknesses (6.12) also agree. 
Jungclaus gives a solution for diverging flow which is equivalent to the case sketched 
in Fig. 3(b); he notes that this solution breaks down for Ry S < 2, and infers that 
separation occurs at this point. However, this result is valid only for solutions of 
boundary layer type, and our solutions (derived from the exact equations) are more 
significant in these circumstances; that is, separation can occur only when Ry S < 3. 
For intermediate values of Ry S, the solution appears to be similar to that given 
by Jungelaus for a jet in a source flow. 

In the discussion following the paper several comments were made concerning the 
validity of the assumption of a small but finite magnetic Reynolds number. It is 
emphasized that in this problem the Reynolds numbers have definite values (namely, 
Q/n and Q/v) in contrast to say, Stokes flow where the appropriate length scales are 
different for flow near the body and at great distances from it. Thus the neglect 
of terms of order R%, in this work probably does not lead to the difficulties found 
in the Stokes and Oseen problems. 

Jungclaus did not show that his solutions can be related to a possible three- 
dimensional experiment; however, our treatment of this aspect of the problem 
should answer any objections that can be raised on these grounds. 

+ G. Jungclaus, Rev. Mod. Phys. 32 (1960) 823. 
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SUMMARY 
The equations of motion for the flow of a non-Newtonian fluid between two infinite 
plates, one of which is rotating and the other is at rest, have been solved approxi- 
mately. It is found that under certain conditions depending on the distance d 
between the plates and the angular velocity Q, the non-rotating plate experiences 
a suction, but when d is decreased sufficiently it experiences a thrust. In the case 
of Newtonian fluids the non-rotating plate always experiences suction. 


1. Introduction 

STEWARTSON (1) has considered the flow of a Newtonian fluid between 
two infinite disks when (a) the disks rotate in the same sense, (b) the disks 
rotate in the opposite sense, (c) one disk rotates and the other is at rest. 
Taking the plates to coincide with the planes z = 0, z = d and the velocity 
in the direction of 6 to be v = rG(z), where (r,6,z) is a set of cylindrical 
polar coordinates, he found in case (c) that for small Reynolds number R 
(= Qd*/v,, where v, is kinematic coefficient of viscosity) G(z) varies 
directly as the distance from the rotating disk and that the series for G(z) 
is rapidly convergent for R < 10. He observed that the boundary layer 
theory cannot be applied for R < 40, for the boundary layer begins to 
grow at the rotating disk when R is about 10. He carried out numerical 
analysis for the case R = 64, when the distance between the disks is 
8mm. In this case the boundary layer is formed near the rotating disk 
and in its vicinity the axial flow is very little affected by the non-rotating 
disk, the modification in the flow near the stationary disk being negligible 
and there being no significant couple on the stationary disk. 

In this paper the equations of motion for the flow of a non-Newtonian 
fluid, confined between two infinite parallel disks, one rotating and the 
other at rest, are solved. The Reynolds number 2 is taken to be less than 
unity and the distance d between the disks is taken to be of order 0-02 mm. 
The solution is expressed in powers of the Reynolds number. It is found 
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that, up to certain values of the distance between the disks and of the 
speed of rotation 2 the non-rotating disk experiences suction but when d 
is decreased sufficiently, it experiences a thrust. In the case when New- 
tonian fluids fill the disks the non-rotating disk always experiences suction. 

The peculiar behaviour of the non-Newtonian fluid can be explained 
by writing the equation of flow in dimensionless form and examining the 
role of different non-dimensional parameters. The equations of motion 
of non-Newtonian fluid governed by the stress and strain-rate relation 
(5) in dimensionless form is 

OW texvd curl w+ |= wait oF ee oF ow 

et LU L?\ da Ox ty fy bz lz 


+ (curl. Vyv-+-4(e0. Vo| = —- 
p 


tials ¢ = Ip|vit— 2 Eel + {er 4 | v |? 
ea} © fey, © jee] 

and v and w are the velocity and vorticity vectors. Hence the flow pattern 
in non-Newtonian fluids depends upon three dimensionless parameters, 
the Euler number EZ = U/,/(2P/p), Reynolds’ number R = LU /v, and 
v,/L*, where L is a typical representative length of the system (in this 
case the distance between the disks). Inspectional analysis reveals that 
cross-viscous effect will dominate the viscous effects if (v,/L?) > R-' or 
K> 1 (K = p,U/p, L). The condition K > 1 can be fulfilled either by 
making U large or L small. For example, when any non-Newtonian fluid 
is sheared between the disks [Popper and Reiner (2)|, one rotating and 
the other at rest, the flow pattern as well as stresses on the disks are 
completely different from those of Newtonian fluids when the distance 
between the disks is appreciably reduced. Similarly, when the fluid has 
to pass through fine gaps, it is to be expected that the cross-viscous effect 
will dominate the flow pattern as well as stresses. 


2. Equation of motion 
The equations of steady motion are 


ou ous” 
plu— +w—— - 
or r, 


Cz 


ov Cv uw 
pi u— + w— + = - 
or cz r 


Ow Cw Or’, 
pl(u— +w—]) — — 
| Or ez a 


where u, v, w are velocity components in the directions r, 6, z respectively. 
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The equation of continuity is 
(4) 


The stress matrix (ri) and the rate of strain matrix (e}) in the case of non- 
Newtonian fluids are related as | Braun and Reiner (3)]. 


i x 


ti = —pdi+p,ej+p, eter, (5) 
where y, is the coefficient of viscosity and py, is the coefficient of cross- 
viscosity. 

Consider the motion of a non-Newtonian fluid filling the space between 
two infinite parallel disks, at a distance d apart, one of which (z = 0) is 
rotating with constant angular velocity Q about an axis (r = 0) perpen- 
dicular to its own plane and the other disk (z = d) being at rest. The 
boundary conditions in this problem are 


w=0 atz 
(6) 
7=0, w= at z = d. 
In this problem the distance d between the disks plays an important 
role in determining the nature of the flow and the stresses on the disks, 
so the length is made dimensionless by introducing a new dimensionless 


parameter Mi 


7 = a 
The velocity components and the pressure are taken to be of the form 
u = rQF(n), vp = rQG(n), w = dQH(n), (8) 


9 


_P(n)}. (9) 


? 


Pp fy, $2) — P( n)+ 7 

Substituting these values of velocity components and pressure in the 
equations (1)—(5), we get 

R( F?—G?+ F’H) = F’+aR(F""*—2F F’ —G")—2Bh, (10) 

R(2FG+-G'H) = G’+2aR(F’'G' — FG"), (11) 


2 
R(HH’') P'+H"+oR(8H'H" —4F’ F”) + ral 20R(G' "+ F' F")—P,], 
a 
(12) 


where R = Qd? vy, and a = v,/d? and a prime denotes differentiation with 


respect to ». The equation of continuity gives 
2F+H’ = 0. (13) 


Equating the coefficients of different powers of r on both sides of (12) and 
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integrating the equations thus obtained, we get 
P—P, = 4RH?+2F—14aRF? (14) 
and P, = aR( PF"? +G")-+A, (15) 


where F, is the value of P at » = 0 and J is a constant of integration. 
Using (10) and (15), we get 


R( F?—G?+ F’H) = FP" —aR(F’?+3G'2+2F F’)—2d. (16) 
For small values of R, a regular perturbation scheme for the equations 
(11), (13), and (16) can be developed by expanding F, G, H, A in powers 


of R. We take 
F = f,+f, Rif, R*+-f, R?+... | 


G = got+g, R+9, R*+9, B+... 

H = hth, R+h, Beet | 

A = Ay +A, R+A, R?+A, R3+ 
Substituting F, G, H, and A from (17) in (11), (16), and equating the 
coefficients of different powers of R on both sides of these equations, we 
obtain the following equations: 


o = 2A, Jo = 9, (18) 
fi = fo +390'+ 2fof 0) + (f3—9—S oho) +2A,, 

Gi = (Joho+ a —2a(fodo—So9o): (19) 
fa = 2h of 1+ 6991+ 2fofit2f ofr) +(2fohi— 29091 +S ohartS iho) + 2r2, 
92 = [2fo91+2h/i Go+Goh+93 ho]—2o[(fo9it+F190)—(fo9i—figo)]. (20) 
fs = Af P+ 2f of2+39P + 690 92+ fof et 2hSit Bfefo+ 

+H(Fi+2fofe—Gi—2GoGatS ohet+S iti tS eho) +2Az 
93 = [2(fo92+Si 91 +290) + (Joho +9 hi +92 ho)|— 
—2a[(fog2+Si 91 +S 290)—(fog2thi9itfego)], ete. 
The equation of continuity gives 
2f, +h, =0 (n=0,1,2 


whilst the boundary conditions reduce to 
je=9, gg=1, h=0 
io =9, og =0, hg =0 
and for n = 1, 2 
In =% 9, =9, hy 
i,=%, g,=90, A, =0 
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3. Solution of the equations 
The solution of (18) and (22) satisfying (23) is 


to — 0, Jo — l—y, ho — 0, Ay 
The solution of (19) and (22) satisfying (24) is 


fi = fon? +4 — 
9, = 9 
h, = ion” -$*+d5n° 
and A, = HL- hy 
The solution of (20) and (22) satisfying (24) is 
0, 


1 ee RE er ee ee ee ee ee ee ea 
7007) 7 367 127 357] 307) 3107 “(167 3677 +$n 


The solution of (21) and (22) satisfying (24) is 


, 13823 . 1 7 ~5 100 ge, ag Se oe 10 
Is 3335507 + $3000" 6000") ~~ 36006 b60071° + sie009" — 67507° + ste07° — sidan + 
is 8737 478439 7 Si 76 58.5.) 917.6. £7.16 
+o — }5$h507 + $388507" — 1507" + dto7* — 007° + fadon® — 1889’ + 1é087°), 
0. 


13823 ..2 6877 _ 3 7 ~6) 109 41 8) 87 .9 1 10 
3772000") 1188000 7 18000") 126000") ~~ 50400") ~T 302407) 108007) 


1_4_ ¢y( 158787 92 478439 oe, 


wee 1 rt mig 4 
1188007 189000) ~~ 567000 7? 300") ~~ 150077 


| 53.6 31 .7_i. 5 8 5 9 
$00” — 9007 + 30a” —a550")s 


6877__| 476549 
and As 792000 7 378000°° 


Using (9) and (14), we obtain 
p = —p, 0] (§H?4+2F—140RF?4P,)— (ORF? vRG?2+))|. 


Since p is everywhere positive in the fluid, ®, must be negative; let it be 
equal to —Q*. The expression for 72 is 


a ees @— 7}: 


folk ale 210+ 50 LA, R4-Ay R2-+Ay RO ~)}. 


Neglecting R® and higher powers, we obtain 


Zz 2 2,2(1 
[rz]. d p,2QQ5— 3 pQ r*(io—%). 
The average normal force on the circular portion of radius a of the non- 


rotating disk is 
a 


“ | 2nr|72),4 dr $pQ2a?(a— bh) — pw, 295. (25) 
7a* 
0 
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Hence the non-rotating disk experiences suction or thrust according as 


lv, , _ 40», Q3 
d 3Qa? ” 
10v, 40v, Q2 


or pet om 1 > Se, 


d* 3Qa* 


The results for the Newtonian fluids can be obtained by putting »v, 
in the above analysis. 

In the Newtonian case the average normal force on a circular portion 
of radius a of the non-rotating disk is 

— (Hyp 2Q5 + sopQ?a*), 
which is always negative. Hence the non-rotating disk will always ex- 
perience suction. 

Reiner (4) performed an experiment in which a gap of air was main- 
tained between two disks, one rotating and the other at rest, and found 
that the manometer placed at the centre of the non-rotating disk indicates 
a thrust. He showed that the above phenomenon can be explained if air 
is an elastico-viscous material possessing an elastic shear modulus. Taylor 
and Saffman (5) explained the above effect on the basis of compressibility 
of air and due to a small error in the perpendicularity of either of the plates 
to the axis of rotation (in Reiner’s experiment). Their point of view is that 
the air is governed by the ordinary Newtonian law of compressible fluid. 
The above phenomenon can also be explained if we attribute the property 
of cross-viscosity to the air and apply the above theoretical analysis. 


4. Conclusion 

When a non-Newtonian fluid is sheared between two disks, one rotating 
and the other at rest, the stationary disk experiences suction when the 
space between the disks exceeds a certain critical value, but the suction 
is changed into thrust when the space between them is less than this 
value. In the case when the space between the disks is filled with New- 
tonian fluid the stationary disk always experiences suction. 
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SUMMARY 

Taylor instability at the interface between two semi-infinite viscous fluids is 

considered. The exponential rate of growth of perturbations is computed for various 
density ratios in the limit as the viscosity of one fluid tends to zero, and is compared 
with the solution of Chandrasekhar (1) for equal kinematic viscosities. 
Tue effect of viscosity on Taylor instability at the interface between two 
semi-infinite fluids was studied by Chandrasekhar (1) assuming the kine- 
matic viscosities of the two fluids to be equal. Tchen (2) made a similar 
assumption, replacing both viscosities by a weighted mean, and approxi- 
mated to the roots of the characteristic equation for the exponential rate 
of growth or decay of perturbations. Instability is most important at large 
density differences, when the viscosity difference is usually large, so that 
these results might not give sufficiently accurate estimates of modes of 
instability. The computations of Bellman and Pennington (3) refer to 
a viscous fluid bounded by a vacuum. 

In the present note the characteristic equation is found for arbitrary 
viscosities, Chandrasekhar’s equation being a particular case. Numerical 
solutions for the exponential mode of instability are then obtained in the 
particular case in which the viscosity of one fluid vanishes. Two cases 
arise, the viscous fluid being either the lighter or heavier of the two. For 
instability, the acceleration, modified by a surface tension term, is taken 
as acting in both cases from the lighter to the heavier fluid. When the 
viscous fluid is the heavier the rates of growth of perturbations are very 
close to Chandrasekhar’s values, there being a mode of maximum insta- 
bility at all density ratios, that is, a perturbation wavelength whose rate 
of amplification is greater than that of all other wavelengths. When the 
viscous fluid is the lighter of the two, modes of maximum instability occur 
only at small density differences, while for large density differences the 
rate of growth of perturbations is close to the inviscid solution of Stokes 
(4) and Taylor (5), as we would expect. A rough estimate of the rate of 
growth of perturbations for any viscosity and density ratios can be 
obtained by interpolation between the present results and those of 
Chandrasekhar. 

Two semi-infinite viscous fluids in contact at the plane y = 0 with 
surface tension 7’ are subject to an acceleration g in the y-direction. In 
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Fic. 1. Rates of growth of perturbations at a viscous interface when 
v = 0 (full curves) and v = v, (dotted curves) for p, > p. 


TABLE 1 


Modes of maximum instability for v = 0 when p, > p (X»,¥n), 
and vy = v, (X,,, Y,,) 
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the region y > 0, the density, viscosity, and kinematic viscosity of the 
fluid are denoted by p, u, v, while in the region y < 0 these quantities are 
denoted by p;, 4;. ¥;. The interface is modulated by a two-dimensional 
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perturbation of the form y = A exp(ikx-+-ot). If the components of velocity 
in the region y > 0 are derived from a stream function 

yb = x(y)exp(ika+-ot), 


then the elimination of the pressure from the linearized equations of 
motion leads to the usual Orr-Sommerfeld equation 


a(x” —k*y) = (x —2k*y"+ ky) (1) 

















Fic. 2. Rates of growth of perturbations at a viscous interface 


when vy 0, for py < p. 


with a similar equation in y, and y, in the fluid occupying the region 
y < 0. The appropriate solutions are 


x = ae~*v + be-*yy, X1 = cekY+-dekny, 


where y? = 1+o/k*y, vi = l+o/k*y,. (2) 


The boundary conditions at the interface are the continuity of both 
components of velocity and both components of stress. Elimination of 
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the constants a, b, c, d from these conditions gives the characteristic 
equation 
Ze y— 1) oy 1) + ply — 1) + ey) + 
+iely+ 1) + pala t Del? + 1+ ealvit+ 1+ (ei—p)g'/ko} = 0, (3) 
where p = pv, uy, = p,v, andg’ = g+Tk*/(p,—p). Withy = v, and T = 0, 
(3) reduces to the corresponding equation of Chandrasekhar’s paper. 
When v = 0, (3) becomes 
2k4p, vi(y1— 1)?+-2k*p, v, o+-0%(p, +p) +-kg'(pi—p) = 9. (4) 
instability occurring when g’(p—p,) > 0. 
In the first case, p, >p and for instability g’ <0. In Fig. 1, 
Y = o(v/g’*)' is plotted against X = k(v*/\g’|)! for various values of 
K, = (p,;—p)/(p, +p). Dotted curves show Chandrasekhar’s values for 
v=v,. The values of X and Y for the modes of maximum instability 
with v = 0 are given in Table 1, and denoted by X,,,, Y, 


m? 


while the corre- 
sponding values for v = v, are tabulated for comparison and denoted by 
ne Bie 

For the second case of instability, p, <p and g’ > 0, Fig. 2 shows 
Y(X) for various values of K, = (p—p,)/(p+p,). For small values of Ky, 
that is for small density differences, Y(X) is again similar to the solution 


for v = v,, exhibiting modes of maximum instability, while for large K,, 


Y(X) resembles the Stokes-Taylor inviscid solution. 
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ON THE OPTIMUM CONTROL OF SHIPS’ STABILIZERS 
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{Received 14 November 1960] 


SUMMARY 

It is becoming an increasingly common practice to fit one or more pairs of fins 
to the sides of ships so as to reduce rolling in a rough sea. The anti-rolling moment 
generated by the fins depends on their attitude, and in all existing systems this 
attitude is controlled, not by a direct recording of the state of the sea, but by a con- 
tinuous measurement of the angle of roll and of the rate-of-roll of the ship. In the 
literature one finds references to the relative amounts of feedback from the roll and 
from the rate-of-roll which have been adopted in particular commercial control 
systems; it appears that these relative amounts have been determined, and may 
be adjusted during the voyage, on empirical grounds, but that no theoretical pre- 
dictions have been made of the optimum combination of these two types of control. 

The special feature of what may appear at first sight as a simple application of 
existing control theory lies in the fact that there is a definite limit, imposed by the 
mechanical design and by the hydrodynamic characteristics, to the deflexion of the 
fins; in other words, the magnitude of the feedback from the ship’s rolling motion 
cannot exceed a certain prescribed amount. On the other hand, one presumably 
desires the feedback to attain its maximum magnitude as often as possible so as to 
achieve the greatest degree of stabilization, and it is the aim of a theoretical in- 
vestigation to determine the feedback coefficients of the roll and rate-of-roll to 
satisfy this requirement. 

So posed, the problem seems intractable for a random sea—that is, with a random 
forcing function in the governing differential equation for the roll angle. But for 
the special case of a sinusoidal forcing function, the solution is simple and is given 
in this paper. The final result, once it is obtained, seems intuitively obvious: the 
roll can be reduced by exactly, but by no more than, the roll angle which the 
stabilizing fins could produce (by a sinusoidal movement over their whole range) in 
a calm sea. The values of the feedback coefficients are given for this optimum control, 
and suggestions are made as to the procedure which might be adopted in a realistic, 


random sea. 


Nomenclature 

time. 
frequency of undamped unforced rolling. 
coefficient of natural rate-of-roll damping. 
angle of roll. 
frequency of quasi-steady motion. 
maximum wave slope of quasi-steady sea. 
maximum equivalent wave slope of fins. 

FW fin effectiveness, assumed less than unity. 
coefficient of roll feedback. 
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coefficient of rate-of-roll feedback. 
amplitude of roll. 

value of R with no feedback. 

value of R with rate-of-roll feedback only. 
optimum value of R. 

Ropt/ Ro- 


~1\2 
nee ane Sd 
b+") ] 
defined in equation (11). 
m defined in equation (19). 


1. Introduction 

Tue design of an arrangement of stabilizing fins for the reduction of roll 
in ships poses two major theoretical problems. First, the fins themselves 
must be shaped so that they exert the required hydrodynamical forces 
with the least loss of energy. This problem is closely analogous to the 
design of aerodynamic lifting surfaces; there is a great mass of aeronautical 
literature relevant to this, and it is not considered here. Secondly, the 
control system governing the attitude of the fins, and hence the forces 
generated by them, must be such as will make the most effective use of 
those forces. 

This paper considers the latter problem, within the limitations of the 
following assumptions. It is first assumed that the angle of roll is governed 
by a single second-order linear differential equation: thus we presuppose 
(a) no coupling between the rolling and other motions of the ship, and 
(6) righting and resistance moments which are directly proportional to the 
angle of roll and to the rate-of-roll respectively. These assumptions are 
those commonly made for small angles of roll, and presumably become 
increasingly valid as the roll is reduced by stabilization. It is also assumed 
that the forcing function—namely, the moment applied to the ship by 
the action of the waves—is sinusoidal. To understand better the extent 
to which results based on this assumption might be applicable to a sea 
of random wave form, we shall briefly interpret the well-known charac- 
teristics of solutions of second-order linear equations in terms of the rolling 
motion of a ship. 

First, consider a ship rolling freely in a calm sea, having initially been 
heeled over by some artificial means. In the absence of any hydro- 
dynamical frictional forces, the ship oscillates in roll with its natural 
frequency f,, and with undiminished vigour, as shown in curve (i) in Fig. 1. 
In reality, frictional moments (which are assumed to be proportional to 
the rate-of-roll) cause a decay in the rolling, according to curve (ii), in 
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practice, the angle of roll decays only slowly which means that the damp- 
ing coefficient ¢ is small and also that the frequency of the damped rolling 
hardly differs from the natural frequency f,. If the curve (ii) were to be 
obtained experimentally, then it is possible to deduce (as shown in section 
3) the values of the coefficients w and € which appear on the left-hand side 
of the roll-equation (1). 





() Caim-sea roiling with no dcmping 
(i) ——— Calm-sea rolling with damping 


Rolling between large isolated waves, spaced at large 
random intervais,with small random waves superimposed 


(iv) initial disturbance of large isolated waves 


Fic. 1. Sketches which illustrate certain conditions of rolling. 


The effect of a randomly distributed wave system may be considered 
by thinking first of random waves of small magnitude among which appear, 
at large random intervals, much larger single waves. A large wave initiates 
the type of motion already sketched in curve (ii) and this is slightly 
disturbed by the small waves; thus between successive impacts of large 
waves, the roll would be of the kind shown in curve (iii). The predominant 
frequency of the ship’s rolling is, even in this case, still roughly equal to 
the natural frequency, though the disturbance produced by each large 
wave as it appears may introduce a transient motion of different frequency, 
as indicated in curve (iv). 

If now we imagine the large waves to form a regular pattern—and in 
particular to act in a sinusoidal fashion—with a frequency which is com- 
parable to the ship’s natural frequency—then the rolling motion takes on 
an entirely different character. The ship now tends to roll with the same 
frequency as that of the large waves, and to abandon altogether its own 
natural frequency; this is roughly illustrated in Fig. 2 in which curve 
(i) shows the sinusoidal large waves with small random waves superim- 
posed, and curve (ii) indicates the ship’s response. 
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While there is, inevitably, a certain artificiality in this simplified dis- 
cussion, we can deduce in a rough way from Fig. 2 that a ship tends to 
respond mainly to some dominant frequency in the wave pattern (provided 
this is of the same order of magnitude as the ship’s natural frequency). 
What is more important is that stabilization will probably enhance this 
effect, for stabilization will result in a much more rapid decay of rolling 
(at the natural frequency) after a large disturbance. Thus in curves (ii) 
and (iii) of Fig. 1, the effect of an isolated large wave will be greatly 
*h 


() 





. 


(i) Sea 
(i) Ship 





Fic, 2. A sketch to illustrate that in a sea with a predominant frequency 
a ship will roll (with a lag in phase) with an equal predominant frequency. 


diminished, and the ship will have a lesser tendency to respond at its own 
natural frequency. On the other hand, for the regular wave pattern of Fig. 
2, although stabilization may reduce the magnitude of the rolling, it cannot 
affect the frequency of either the wave pattern or of the ship’s response 
to it. Thus in a stabilized state, a ship is more, rather than less, likely 
to respond to some dominant frequency in the sea. Of course, for a truly 
random sea, it is impossible to discuss matters at all on such general lines. 

The analysis of this paper concentrates attention on an idealization of 
the state of affairs in Fig. 2; that is, it assumes (a) a regular wave pattern 
which produces a sinusoidal rolling moment on the ship, and (6) that the 
ship has reacted to this wave pattern for a time sufficiently long for the 
effects of all initial disturbances, such as those shown in Fig. 1 (iv), to 
have been damped out. Under these circumstances, both the ship and 
the stabilizing fins will move sinusoidally at the same frequency as that 
of the waves, though in general all three movements will be out of phase 
with each other; and motion of this kind will be referred to as ‘quasi- 
steady’. 

We referred earlier to the two types of feedback which have been used 
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in fin-control systems, namely, feedback proportional to the roll angle and 
to the rate-of-roll. It is possible that advantages could be obtained in real 
seas by the empirical use of roll-acceleration feedback; we ignore this 
possibility in the subsequent analysis since, in the particular case con- 
sidered, the roll acceleration is directly proportional to the roll-angle and 
so acceleration feedback can contribute nothing of a different kind. 

To sum up the preceding discussion, we would hope to justify the follow- 
ing analysis of the quasi-steady rolling motion by remarking that most 
wave systems in deep water have a dominant component of frequency, 
and that in a highly-stabilized motion, the ship will have a strong tendency 
to take up that same frequency. 


2. The equation of roll and the general quasi-steady solution 
Let the natural frequency of the ship’s rolling be w/27, the frequency 
of the sinusoidal wave system be ow/2z7, and the maximum slope of the 
wave be W. Then the equation for ¢, the angle of roll, may be written 
in the form 
1 d*p | dd | 


wd?" w dt 


$ = W(14 2oif)exp(iwot) kee 
a 


in which ¢ is the coefficient of roll-damping, and K, and K, are the feed- 
back coefficients of roll and rate-of-roll respectively. This form of the 
equation follows that given by Chadwick (1) except for the introduction 
of the complex wave-forcing function W(1+-20if)exp(iwot) which is more 
convenient than its real part. 
The maximum angle of roll R may easily be found by the substitution 
¢ Rexpliwo(t x)] (2) 
into the roll equation (1); here, « is the phase-lag between the wave 
oscillation and the ship’s roll. By taking the modulus of both sides of the 
resulting equations, we find 
R| —o*+-io(2¢+ K,)+ 1-4 k, W \1+-2t0% |, 
1+402f2 
(1—o?+ K?)?+-07(K,+ 2¢)? 
which determines the roll angle R in terms of W, a, (, K,, and K,. Of 
these five quantities, 


or R? Ww? (3) 


(a) W can be regarded as fixed; or otherwise, one can think of the ratio 
(R/ W) which, from (3), expresses the magnitude of the roll; 


(6) Cis a characteristic quantity for a ship; its experimental determina- 
tion is discussed briefly in the next section; 


(c) o is a parameter describing the frequency of the wave system in 
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terms of w; w/27 is the natural frequency of the ship’s rolling which, like 
{, can be determined experimentally; 


(d) K, and K, are the feedback coefficients which it is our object to 
choose so that R is minimized. But K, and K, are also conditioned by 
the fact that the maximum feedback is pre-determined (as was discussed 
in the summary) and the implications of this are further discussed in 
section 4. 


3. The experimental determination of w and ¢ 

These two quantities may be determined by releasing a ship from an 
artificial list in a calm sea, and measuring its subsequent rolling motion 
to produce a record such as that in Fig. 1(ii). In the absence of waves 
or feedback—that is, when W = K, = K, = 0—the solution of equation 
(1) is 

¢ = exp(—lwt)[ Bexp{—iwt,/(1—¢*)}+Cexpfiwt/(1—*)}], (4) 
where B and C are constants. Thus the amplitude of roll is decreased, in 
each complete oscillation, by the factor 5 given by 
5 = exp[—2z7f(1—{)-*]. 

5 is usually of the order of 0-7 and so ¢ is of the order of 0-07; thus a good 
approximation to 6 is exp(—27¢) and so 


] i e 
(= = log (;). (5) 


Also, the frequency of free rolling may be taken as w/27, from equation 
(4). Values of w and { may thus be determined from free-rolling trials. 

The value of { so obtained may be checked with that given from 
observations of the ship’s motion in a wave system whose frequency equals 
the ship’s natural frequency w/27. This is called the synchronous state 
in which rolling of large amplitude may be expected. In this case, we 
have o = l and KA, = K, = 0, so that from equation (3), 

R/W = [(1+(1/20)7}*; 

if R/W can be estimated experimentally, then a value for { is deduced 


which may be compared with that derived from calm-sea trials, 
In what follows, w and ¢ are regarded as known quantities. 


4. The fin effectiveness E 

We have previously remarked that the attitude of the fins is limited; 
on the other hand, it is clearly advantageous to use the fins to their fullest 
extent. Now in the quasi-steady motion, the fins move in a sinusoidal 
fashion and it is desired that they move in this way over their full range. 
How is this requirement best expressed ? 





ON THE OPTIMUM CONTROL OF SHIPS’ STABILIZERS 369 


It is best to consider this question by reference to the method by which 
the effectiveness of the fins might be measured in practice. By analogy 
with the procedure for determining w and C, the fins could be set to exert 
their greatest rolling moment and the steady angle of roll then measured. 
Let this angle be F. Then F will be a characteristic angle for the fin 
system which can be determined for any speed of the ship. It corresponds 
to the angle W, the maximum wave angle, which, if maintained steadily, 
would result in a constant roll-angle of W. Thus we propose to measure 
the fin effectiveness EZ in any given sea by the ratio of F to W. 

Now from equation (1) we see that if the fins in a calm sea will produce 
a roll equivalent to a wave system of slope F,, then 

|K, dd 


F(1+ 20i1f)exp(icvot) = | +K,¢. (6) 
lw dt 


Substituting the expression (2) for ¢, we have 
F*(14-40°f?) = R%(o®K?+ K?), (7) 
and substitution for R? from equation (3) gives 
E*{(1—o? + K,)*-+-0%(Ky+2L)?] = (K?+02K3) (8) 
where E=F/W. (9) 


Equation (8) is therefore to be interpreted as the functional relationship 


between K, and K, to ensure that the fins are used to their full capacity. 
E is a parameter in this relationship and equals the known characteristic 
F of the fin system at the speed of the ship under consideration divided 
by the maximum wave slope of the sea. We shall find later, in section 7, 
that there is no purpose in allowing KE = F/W to be greater than unity, 
and so we shall in fact assume henceforth that 


E= FiW <1. (10) 


E is called the fin effectiveness. 


5. Rate-of-roll control only 

Before finding the optimum control, it is interesting to establish the 
effect of a rate-of-roll control alone, since this has been used widely in 
practice; and we shall assume, optimistically, that this control is adjusted 
for the best results, namely that the fins are used to their full extent accord- 
ing to the condition (8). 


Let us put, therefore, 
K, . K, = 20A. (11) 
Then, from equation (3), 
‘ 1+ 407? 
R=- B= W* —_._.___- — 
” (1—o?)? + 407f7(1 +-A)? 
Bb 


(12) 
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where 2, is the maximum roll under this control, while from equation (8), 
E%(1—o)?-+-40°f%(1-A)*] = 40°f2)2. (13) 


To obtain the roll-angle R, directly in terms of EZ, A is to be eliminated 
from equations (12) and (13); after some algebra, we find 


i (Ste) 





_ Se 2 al 


where Pie [p+ (*5¢-)} j 


Tr 


(14) 


This formula therefore gives explicitly the angle of roll under rate-control 
only, if the fins are used to their full effectiveness as given by E. The roll 
so obtained is compared in section 7 with the optimum roll to be discussed 
in the following section. 

For no stabilization, the fin effectiveness Z is zero, and equation (14) then 
reduces to R R, 


— So (I +40°%?)t 2 ( 15) 
W W — {(l—o*)?+40?}! 
This result also follows from equation (3) with K, = K, = 0. It may easily 
be verified that, as the fin effectiveness increases, the roll R,, as given by 
equation (14), decreases; while for E = 1, R, apparently becomes zero. 
This, at first sight, appears a reasonable result; for the fact that E = | 
means that the fins alone can produce a steady roll exactly equal to the 
maximum slope of the waves, and so it should be possible to set the one 
against the other all the time and thus obtain zero roll. But it may easily 
be verified that this state of affairs involves an infinite feedback coefficient 
K, and so is unrealistic. In different words, one cannot use an identically- 
zero quantity to be multiplied as a feedback signal. 


6. The optimum combination of roll and rate-of-roll feedback 
The condition that the roll R given in equation (3) is a minimum is 


d{(1 —o?+K,)?+0%(K,+2f)2] = 2(1—o®+ K,) dK,+20%(K,+2L) dK, = 0. 
(16) 
This must be coupled with the condition (8) which in differential form is 
E*d[(1—o? + K,)?+-07(K,+ 2f)?] = 2(K, dK,+0°K, dK,). (17) 
In order that (16) and (17) may be simultaneously satisfied, we require that 
20K, = (1—o?)K, (18) 

which is obtained by eliminating dK, and dK,. 


Equation (18) specifies the proportion of roll and rate-of-roll feedback 
to the fins. This proportion (which may be regarded alternatively as 
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determining the phase difference between the roll and the fin movement) 
depends only on o, that is, on the ratio of the frequency of the motion 
to the ship’s natural frequency; it does not depend on the fin effectiveness 
or on the speed of the ship (except in so far as the speed determines the 
frequency of wave encounter). 

In the case of synchronism, o = 1, and we see that K, = 0; in other 
words, the pure rate-of-roll control is the optimum in this case, with the 
fins moving exactly 90° out-of-phase with the waves. However, the wave 
frequency only has to depart by more than about 6 per cent from the 
synchronous frequency for K, to be greater than A,, that is, for the roll 
control to be of the greater importance. 

The actual values of K, and K,, rather than their ratio, are found from 
equations (8) and (18). From (18), putting 

a, 3 BL, (19) 


equation (8) gives 


_ BY (1 —o?)? +-40°C?]( 1 +a)? = [(1 —0?)? + 40°C? ]p?. 
E 
+i—-z 


Now in choosing the correct sign, we note first that in equation (1) K, must 


Hence ph (20) 


be positive (or at worst greater than —2{); otherwise there is the possi- 
bility of unstable oscillations. Thus, from (19), ~ > 0 and so, from (20) 
not only must the positive sign be taken (since £ is essentially positive) 
but also £ must be less than unity. This at first sight seems a surprising 
requirement since it means that the fin system must always operate at an 
equivalent wave-slope less than that of the wave system it is intended to 
uateract. The paradox can be resolved as follows. 
First, the values of K, and K, are, from equations (19) and (20), 
E » 
"7a alld a“), kK, i 
The minimum roll-amplitude is then given from equation (3) as 


1 + 4072? Some » 
R anne | (1 — #). 
opt (1 aan ( 


kK, 


The two terms on the right-hand side implied by the binomial (1 — £) are 
precisely the roll-amplitudes generated first by the sea alone, and second 
by the fins alone. Thus the physical interpretation of the vuptimum control 
is that the fins reduce the roll by exactly the roll which they themselves 
could produce in a calm sea. This is, of course, a consequence of the correct 
phase relationship implied in (18); any other phase difference between ship 
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and fins would result in a greater roll. The result may be illuminated 
further by the remark that, with the optimum vaiues of K, and K, given 
in (21), the feedback [(K,/w)(d¢/dt)+K,¢] may be shown to be exactly 
in phase with the waves; thus in the optimum control, the fins are moving 
exactly in step with the sea. 

Returning now to the restriction Z < 1, we see that the residual roll 
R,,, can be as small as we please for this range of values of Z. Of course, 
for values of E close to unity, the feedback coefficients given in (21) become 
large, and practical difficulties will arise in the control system; and the 
limiting case HE = 1—of infinite feedback and zero roll—can have no 
practical meaning, as has already been discussed in section 6. Neverthe- 
less, the important result remains that, with a sufficiently sophisticated 
control system, the roll can be reduced to any level for the range of fin 
effectiveness given by EH < 1. 


7. Some general remarks 
It is interesting to compare the optimum roll with that due to the rate- 
of-roll control only. From equations (14) and (22), we obtain at once 
pa Bo IB) BH LY 
R, -E+[H*+77(1— E?)}* (1+ 2) 


10 


(23) 
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Fic. 3. Curves showing the ratio of the rolls in the optimum and in the 

rate-control cases, as a function of the frequency parameter 7, for various 
values of the fin effectiveness E. 


This ratio is plotted in Fig. 3 for 1 < + < 7 and for a number of values 
of E. At + = 1, that is for synchronous rolling with o = 1, the optimum 
control is identical with the pure rate-of-roll control, as was mentioned in 
the previous section; but as o departs from unity (or as 7 increases from 
unity), the optimum control becomes steadily more advantageous. For 
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either large or small values of o (compared with unity), r takes the 
asymptotic value [(1— #)/(1+ £)}}. 

The graph confirms that, for the greatest advantage, EZ should be as 
near unity as practicable. Thus, as a practical requirement, the fins should 
be able to maintain a keel angle nearly equal to the slope of the largest 
waves likely to be encountered; it is at this stage that the hydrodynamic 
design of the fins has to be considered, though it is outside the scope of 
this paper. 

The analysis of section 6 is not susceptible to extension to the case of 
a random wave system. For the reality of a sea, a dominant frequency 
must be obtained so that the ratio of the feedback coefficients may be 
evaluated from equation (18). As an operational procedure, pure rate-of- 
roll control could first be applied; as was argued in the introduction, the 
ship would then, in this damped condition, roll with a frequency associated 
with the sea, rather than at its own natural frequency. With a value of 
o so determined, the optimum feedback ratios could be set, the actual 
magnitudes being adjusted to ensure full use of the fins. Thus in practice, 
it is doubtful whether use would be made of the results (21); the necessity 
of assessing a value of W for a random sea, and hence of #, would hardly 
arise. 

Finally, we may mention that there is no difficulty in incorporating a 
time-lag 7’ of the servo-mechanisms into the analysis, provided that it is 
a constant. The effect is to insert a factor exp(—iwoT’) outside the square 
bracket in equation (1). The final result, then, for the values of K, and 
K, is E 

K, —.. {cos wa T' —[20f/(1—o*) sin woT'}( 1 —o?) 


1—E 


kK, 


{cos wa T' +-[(1—o*)/20€]sin woT'} 22 


E 
1—E 


which may be compared with the values given in (21). 
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SUMMARY 


In the first part of the paper, the temperature at all points of a semi-infinite solid 
heated at a rate Q(y,t) over a strip —h < y < A of its surface z = 0 is obtained 
from the analysis of Lowan (1). It is shown that with a constant heating rate, the 
ratio of the mean temperature at the surface strip to that when the same heating 
rate is applied to the whole surface —« < y < © is 1—(1/2h),/(xt/zr) (« = thermal 
diffusivity, ¢ = time), provided the diffusion depth 2,/(«t) is moderately small. This 
permits one to estimate the error introduced by edge effects in the calibration by 
short-duration Joule heating of thin film resistance thermometers for use in shock 
tubes. 

The remainder of the paper deals with the variation of the surface temperature 
T(t) at the stagnation point of a bluff body in a shock tube. It is noted that the 
heating rate Q(t) may be written 


(pho tQiey = {nO BE) (0 <<} 


\A(T,—T(t)) (t > to), 

where A,, A,, 7, T, and (pke)~* are constants, and ¢, is the time that elapses between 
the passage of the shock wave and that of the interface between test and driving 
gases. The exact solution for T(t) in these circumstances is found by operational 
methods in section 4, and expanded to give a simple approximation for smail values 
of (t—t,). 


1. Introduction 


Tus note contains the solutions to two unsteady heat-conduction prob- 
lems that have applications in experimental work in shock tubes. 

The first is a two-dimensional problem: heat is supplied at a rate Q(y, t) 
to an infinite strip —h < y <h in the surface z = 0 of a semi-infinite 
solid, and an expression is obtained for the temperature 7'(y,z,?) at all 
points of the solid. The method for finding this was given by Lowan (1), 
but there is an error in the later stages of his solution. In the same paper, 
however, he gave the correct solution of the more general problem in 
which Q is a function of x, y, t. This is quoted in section 2, and used in 
section 3 to deduce the correct result for the two-dimensional case (which 
Lowan himself treated ab initio). 

From this, the mean temperature in the surface strip is obtained as a 
function of time, and it is shown that with a constant heating rate, the 

(Quart. Journ. Mech, and Applied Math., Vol. XIV, Pt. 3, 1961) 





376 D. A. SPENCE 


ratio of the mean surface temperature when heat is applied to the strip 
only, to that when heating is uniform over the whole surface z = 0, is 
1—(1/2h),/(«t/r), where 2,/(«t) is the penetration depth at time ¢. This 
has an application in the calibration of thin metallic resistance thermo. 
meters for measuring short duration aerodynamic heat transfer rates. 
Such gauges, which over short times approximate to semi-infinite strips 
(typically, the length is ten times the width) are used under one-dimen- 
sional heat flow conditions. They are calibrated by Joule heating, which 
does not reach the retnainder of the surface (except possibly by radiation 
and free convection, which are thought to be very small effects). The con- 
ditions are therefore approximately two-dimensional, and the result now 
obtained permits one to correct the calibration for edge effects. 

The second problem arises in connexion with measurements of the 
temperature at the stagnation point on a blunt body mounted in a shock 
tube. For the heating rate at such a point (in undissociated gas) Fay and 
Riddell (2) give the formula 


Q = 0-94A(p,, 14)" (Po Ho)” (ho—h,,).- 
(The presence of dissociation would introduce a further factor, involving 
the Lewis number.) Here, ho, pp, 4, are the enthalpy, density, and viscosity 
at an isentropic stagnation point; /,,, p,, 4, those at wall temperature; 
and A? is the velocity gradient at the stagnation point. Besides giving 


heating rates under test conditions for comparison with this formula, the 
temperature record provides a means of identifying the interface between 
the test and driving gases, and of recognizing combustion when it occurs 
in the driver. In a typical case, the test gas behind the initial shock wave 
might be at 6,000° K, and the surface temperature of a model could be 
raised to 1,000° K during its passage, which occupies several hundred 
micro-seconds. The test gas is normally followed by cooler driving gas, 
whose temperature of perhaps 1,000° K may be comparable with that of 
the surface; another possibility is that combustion in the driver may have 
raised its temperature to something like that of the test gas. 

In either case, it is of interest to calculate how the surface temperature 
should vary if the heating rate changes abruptly after a certain time fy. 
The variation with wall temperature of the term (p,,,,,)°! in the Fay and 
Riddell formula is trifling, and for wall temperatures 7’ up to about 
1,000° K, the enthalpy h,, may be written as C7, with C, taking an 
appropriate constant mean value over the range. The formula is thus of 


the form (pke)-*Q(t) = A[7,—T (0). 


(The density p, thermal conductivity k, and specific heat of the model c, 
have been introduced on the left-hand side for convenience.) We investigate 
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in section 4 the case where this holds for t < t,, after which A, and 7, are 
replaced by different constants A, and 7,. In practice, the test gas is so 
hot in the initial stage that 7, > 7(t) and the temperature variation of 
the surface could be ignored, but it is possible to solve the full problem 
specified here without much extra complication. The same methods would 
also permit 7, to be a prescribed function of time—as is in fact the case 
when the gas properties vary as a result of shock attenuation—but this 
has not been done. 

The solution found in section 4 can be expanded to show that for small 
{—t,, i.e. just after the passage of the interface, 


T)—Tr) = 210.9. —2}* + 00-t), 
mpke 


where Q, and Q_ are the heating rates immediately after and immediately 
before f). This means that in general the temperature will commence to 
fall parabolically with time after the passage of the interface, unless com- 
bustion is taking place in the driver when a parabolic rise would be 
expected. In practice the interface is not completely abrupt, being smeared 
out to some extent by diffusion of both heat and mass. A recent investiga- 
tion by Goldsworthy (3) shows, however, that these effects are confined 
to a narrow contact layer of width ~ ,(«,t), where «x, is the thermal 
diffusivity in the gas, and the qualitative behaviour deduced from the 
assumption of an abrupt discontinuity is likely to be correct, at least for 
short running times. 


2. Temperature in a semi-infinite solid heated at the surface z = 0 
The heat diffusion equation is 
3 
( «vr == U, (1) 
ot 
where 7'(x,y,z,t) is the temperature (measured from its initial constant 
value) and « the diffusivity k/pe (k = thermal conductivity, p = density, 
c = specific heat). At the surface 
(eT 0 t <0) 
(; \ = | : ( (2) 
GzJo.0 \Q(z,y,t) (t > 0) 
is the inward flux of heat per unit area. Using operational methods Lowan 
has obtained the solution of (1), subject to (2), as 
T (x, y,z,t) 


t x x 
] a e o*/4at—r) P [ F (x £)? 1 (y n)* 
inictk | (t—7)h Ir Q(E, , rT) xp| — dnlt—r) déd». 
0 -—-o —-« (3) 
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This is equation (C 9) of his paper. The results for two- and one-dimen- 
sional heat inputs are special cases: 


Two-dimensional heating 
If Q = Qy,t) only, on integration with respect to é 


t oo 
1 e-2*/Ant—7) (y— n)® 
Ty,2.t) = =~ | "ar [ Qn, rhexp| —{ dn. (4 
(y 2.0 2r7k J t—r id 3 An “exp| eal 7 *) 
0 — 2 
Lowan obtained this result ab initio in section D of his paper for the half- 
plane case in which Q = 0 for y < 0. (The factor outside the integral is 
given incorrectly as 1/2k7* in the paper.) 


One-dimensional heating 
If @ = Q(t) only, integration with respect to 7 yields the well-known 
expression , 


] | Q(r)e-2*/4xtt—) 
— dr 


— Gepkey J (t—7)h "7 
0 


(2, 
for the temperature a depth z below a uniformly heated surface. 


3. Heat applied to the strip |y| < h 


Case A of Lowan’s paper is that in which heat is supplied to the strip 
y| <A symmetrically about y = 0, so that 
Qly,t) = ol yt) (Iy| <A) (6) 
(\y| > A). 
The inner integral in the expression (4) for the temperature distribution 
in this case is 


[ Q(n, r)fe-w- Axit—7) 1 9—y+n? axt—7) dn. (7) 
0 
In Lowan’s result, derived ab initio, h takes the place of » in both ex 
ponentials, which is plainly incorrect. It may be noted, however, that in 
the limit A — 0, if 


h 
lim | Q(n, t) dn = Q(t), 
h->0 a 


either expression leads to the known solution 


t 


l e—u*+ 2*)/4x(t—r) d (8 
27k ik i fame Qolr) dr 5) 
0 


for the temperature in a semi-infinite solid heated by a line source on the 
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x-axis. Lowan had satisfied himself of the correctness of his result because 
it reduced to (8) in this limit. 


3.1. Uniform heating over the strip: mean surface temperature 
Setting Q(y,t) = Q(t) in (7), we can integrate with respect to » and 
compute the mean surface temperature as 


h 


_ l 
Teen = ah i Ty, 0, t) dy 


h 
h 
WSS QQ) 2 [ (e RR _ 4 erf — h onl 
(mpke)* J (t —r)i g 2, {x(t—7)} 2,/{«(t—r)}/ 4h 
=h 
Writing h/,/{«(t— és = B, the inner integral is 
B 
3 | 4 (erfu) du = erf B—- 


0 


l—e* hod Baa ie p* 
mB ee: ~ 7p onips 


For 8 > 2 the last term is < 0-00067; since in practice we are concerned 
with cases for which h/,/(«t) > 1, we shall use the first two terms only, 
expressing the mean temperature as 


Tinean = “sais | Oe Nac ) — iF) | 


In particular, if Q is constant 


Tiesm a 29 


ae a. 


Thus the temperature at time ¢ is reduced below that which would be 
observed under one-dimensional heating conditions at the same constant 
rate Q by a factor 1—(1/2z!)(penetration depth)/(width), where the 
‘penetration depth’ at time ¢ is 2,/(«t). The form of this factor was to be 
expected on dimensional grounds; the present analysis has contributed the 
coefficient 1/27'. In practical cases this can only be regarded as an 
estimate, in view of the possibility mentioned earlier of heat reaching the 
surface of the region |y| > A by convection and, in extreme cases, by 
radiation. 


4. One-dimensional heating at rates depending linearly on surface 
temperature 
The remainder of the paper will be concerned with the temperature 
T(t), say, at the surface z = 0 under one-dimensional heating conditions. 
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This is related to the heat input Q(t) by equation (5), with z = 0, and we 
shall find the solution for the case in which 


c _ {A[TZ,—TW)] (0<t <t,), 
k 4 #) = 1*1 
(pke)*Q(t) ya 


the physical significance of which was discussed in section 1. On setting 
z = 0 in (5) we obtain 


(12) 


t 


(mpke)! T(t) = | (t—r)-*Q(z) dr, 
which is written operationally as : 
pT (p) = (pke)*Q(p) (13) 
where 7'(p) denotes the Laplace transform { T(t)e-” dt, and Q(p) is that 
a 


of Q(t). 


(i) For 0 < t < t, the solution is unaffected by the change in Q(t) after 
ty, and may be found by using the value 


x 


(pke)+Q(p) = A, [ [T,— T(t)]e-' dt = d, E — Mp)] (14) 


. 


0 
in (13). Eliminating Q(p), 

T(p) = A,T,/p(p* +), (15) 
whence on inversion we find the well-known solution [Carslaw and Jaeger 
(4) 55), T(t) = T[1—eMerfe(A, vt)] = T1— FA, »t)}, (16) 
say. Clearly, the limiting value 7, is approached at large times, and since 
F(é) = 1—2€/m*+ €-+..., the small-time behaviour is 


\ 


T(t) = 2, 7,( ata (17) 


t 
7 
the leading term being exactly the parabolic temperature rise produced 
by a constant heating rate A, 7;. 


(ii) For t > t), the equation corresponding to (14) is 


te ’ 
(pke)-#Q(p) = | (AT —AT)—Oy—AyT Oem dt-+a,| "2 Tp)], (18) 


0 


the integral from 0 to ¢, being calculated from (16). Using (13) on the left- 
hand side, one finds 


D(p) = (wp+Aq)""[Ae/p){T,—(T, — Tye} + (Ay —Ag)T (p—Aj) *46(p), (19) 
where $(p) = 1— F(A, vte-P* — (A, perf (./(pt,)). (20) 


The first term on the right of (19) is inverted using (15) and (16), together 
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with the fact that #-(f(p)e-”) = f(t—t)). The second inversion is 
carried out in the Appendix, and combining the two results, the tem- 
perature can finally be written as 


m T 1-9 F(A “F 

T(t) = ieT y, 4 »| aonb ‘0)) a I 
0 

A A 


AT, +A, 7, 
ie Ne EE 


ea) Pvt) +( 


yo. 2D 
where 7, = T(t.) = T,[1— F(A, vt)], and J(t) is defined by equation (A 9) 
of the appendix. As too, it is clear physically that the zero heat- 
transfer temperature 7, must be approached, and this is confirmed by 
(21), since F(co) = J(ao) = 0. It may also be verified on setting ¢t = 0 
that the solution reduces to that given by (16) with A,, 7, in place of A,, 7}. 


Expansion for (t—t,) small 

Using the expansion already quoted for F in terms of its first argument, 
together with that of J for small (t—t,) given by (A 10), we find from (21) 
after some algebra that 


T(t)—T, 


0 


dT | 
: + € i (22 
Som ec) + (ar) a —cEp (Q,—@_) to) +-O(t—ty)® (22) 


for small values of t—t,, where Q_, Q, denote the heating rates 
(pke)'A,(T,—T)) (¢ = 1,2) immediately before and immediately after f,, and 


dT l 
Pe . — A, F(A, vt 
(a) saa 0) 


is the derivative of (16) at ft). 


5. Concluding remarks 


From (22) we see that after the passage of the interface the temperature 
will begin to rise or fall parabolically with time according as Q, 2 Q_. 
If the two rates are equal at ft), their difference at later times will affect 
the term of order (t—t,)!, but the gradient d7'/dt will be continuous at t,. 
Fig. | illustrates the kind of variation to be expected from equation 
(22), the constants used in plotting the curves having been chosen 
arbitrarily. 

It may be noted that the result (22) could also be found rather directly 
from an iterative solution of the integral equation for 7(t) near t. We 
may also note that the methods of section 4 could have been applied to 
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T,=1 
A, =Ae* | 
t,= 0.25 


(arbitrary units) 


t, 





i Yy ; 
re) 4 2 ' 3 4 





Fic. 1. Example illustrating effect of a sudden change in heating rate. 


the more general problem in which 4,7, is a known function, say 
(15) then becomes 7'(p) = A(p)/(p'+A,) with the solution 


(nz) —A, F(A, vr) | dr 
V 


t 
T(t) | Ac ‘) 
: 


of which (16) is a special case. 


APPENDIX 
In this appendix the value of 
$(p) 
?—)32) —- ~ 
(Ai—4) (wp +A_)p—Aj) 
will be found for ¢ > ty, d(p) being the function defined by (20). The required 
inverse may be written 
(i a) ( vp—A, 
7 (p Ai)(p 
t 


x 


= <= b(p)e”' dp, (Al 
Ori PP? z ) 
: 

where c is such that all the singularities of the integrand lie in rep < c. It is 
sufficient to take c > 0, since the integrand is bounded at A}, A2. For large values 


of |p| the integrand is O(e?''-%)), so for t > t, the path of integration can be deformed 
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to the lower and upper sides of the negative real axis together with a circle of 
vanishingly small radius described anti-clockwise round the origin. 

Setting p = xe~'* the contribution from the lower half of the axis (described from 
left to right) is found to be 


on Fay ae aa 
Qi | aw 24 yp) (ive + Andie e*! der, (A2) 
0 


and that from the upper half (right to left) is the conjugate of this. The path round 
the origin makes a vanishingly small contribution, since the integrand is bounded 
as |p| > 0. Moreover, ¢(ze~'”) = ¢(xe'”), since 
rte 
. . , 1 fedv 
ierf(—iv(2t,)) ierf(iv(aty)) — | me 
0 
and summing the contributions, the inverse is 
rte 
l A, [ er dv 
me am! f . _F rte) i —atdy 
r+Aj 2 >) [20 P(A, sto)e y qr | ot 1° dz. (A3) 
0 0 
The first part of this integral is found using the standard result [Erdélyi (5) 136] 
| xi(x+A*) et dx 
7 * 
0 
where, as before, F(Avt) = e“erfc(Avt); setting A = A,, A 
from the terms in curly brackets is found to be 


—A, F(A, vt) +A, F(A, vt) + F(A, vtp)[Ay F(A, V(t—ty)) Ag F(Ag V(t—ty))). (A 4) 


It remains to calculate two integrals of the type 


1 


arty -AF(Avt), 


in turn, the contribution 


at, 
- 


rt | et’ dv 
x 


0 
To do this, we note that 
to 
e’ dv 


vd 
on integrating by parts, so that 


.-Foe t 
le-3t /{ 4 } du + constant; 
T. uN \u—ty 
t 


x 


and by letting ¢ tend to infinity, the constant is seen to vanish. Thus the remaining 
terms in (A 3) contribute the term 


u’ 


A " {{ te du 
= | du—t)/(- , (A5) 
t 


where u(y) = e~t¥— ei. (A 6) 


The required integral is the sum of (A4) and (A5). This may be put in a more 
convenient form by noting that 


. . 3 / t du 
F(A, vt)— F(A, vt) =f blu 0) (- ) : (A7) 


u 
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(A 3) is then found as 


— 2A, P(A, vt) + (Ay +Aq) P(g Vt) + P(A, vbg)[Ay F(A, y(t—ty)) Ag F(Ag V(t — ty) + Az F(t), 
(A8) 


x 


Pree ae (/( A _ de 
where Jit) = = {| (_—) : (5) oe \—. (AQ) 


t 
Clearly J(t,) = 0, and the behaviour for small values of (t—t,) is found from the 
derivative 
ad ry [s Zs du 1 F abv) dv A,—Az 
(ah s = | ao \z = pag hfe) u Inne, | yO ae 


te 7) 
integrating by parts. For small values of (¢—t,), therefore, 
J = —(A\—Aq)t—ty)/V(r1ty). 
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(A 3) is then found as 
2A, P(A, vt) + (Ay +Aq) F(Ag vt) + F(A, vtg)[Ay F(A, V(t—ty)) Az F(Ag y(t —ty))]+Aq S(t). 
(A8 

iff / 


du 
where Jit) 


t j lo | ye 9) 
FAP ) J Ne i.) iu * ae 


t 
Clearly J(t,) 
derivative 


x x 
ed eT ie ifs du l * U(v) dv A, —A, 
(ey (4 oe oe | : 
Ct /t=t, 7) ldtn \u—t! jit, 
t, 


u Q7Nlby | vt vit, 
0 


0, and the behaviour for small values of (t—¢,) is found from th+ 


integrating by parts. For small values of (t—t,), therefore, 


J (Ay —Az)(t —ty)//(a1ty). 
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